From 4efe00d8f5a63c5c8cd939eaf588241890fc84a8 Mon Sep 17 00:00:00 2001 From: npriyadarshi Date: Wed, 15 Mar 2017 20:38:50 +0530 Subject: [PATCH 1/4] Corrected the tests in 00_nmatrix_sec.rb and inverse_exact methods --- ext/nmatrix/math.cpp | 80 ++++++++++++++++++++--------------------- spec/00_nmatrix_spec.rb | 6 ++-- 2 files changed, 43 insertions(+), 43 deletions(-) diff --git a/ext/nmatrix/math.cpp b/ext/nmatrix/math.cpp index a69c8f52..880c1cd4 100644 --- a/ext/nmatrix/math.cpp +++ b/ext/nmatrix/math.cpp @@ -188,7 +188,7 @@ extern "C" { // Math Functions // //////////////////// -namespace nm { +namespace nm { namespace math { /* @@ -232,8 +232,8 @@ namespace nm { DType* a = reinterpret_cast(storage->a); IType col_pos = storage->shape[0] + 1; if (M == 2) { - if (ija[2] - ija[0] == 2) { - *result = a[0] * a[1] - a[col_pos] * a[col_pos+1]; + if (ija[2] - ija[0] == 2) { + *result = a[0] * a[1] - a[col_pos] * a[col_pos+1]; } else { *result = a[0] * a[1]; } } else if (M == 3) { @@ -258,7 +258,7 @@ namespace nm { rb_raise(rb_eArgError, "some value in IJA is incorrect!"); } } - *result = + *result = m[0][0] * m[1][1] * m[2][2] + m[0][1] * m[1][2] * m[2][0] + m[0][2] * m[1][0] * m[2][1] - m[0][0] * m[1][2] * m[2][1] - m[0][1] * m[1][0] * m[2][2] - m[0][2] * m[1][1] * m[2][0]; @@ -270,13 +270,13 @@ namespace nm { } /* - * Solve a system of linear equations using forward-substution followed by + * Solve a system of linear equations using forward-substution followed by * back substution from the LU factorization of the matrix of co-efficients. * Replaces x_elements with the result. Works only with non-integer, non-object * data types. * * args - r -> The number of rows of the matrix. - * lu_elements -> Elements of the LU decomposition of the co-efficients + * lu_elements -> Elements of the LU decomposition of the co-efficients * matrix, as a contiguos array. * b_elements -> Elements of the the right hand sides, as a contiguous array. * x_elements -> The array that will contain the results of the computation. @@ -291,7 +291,7 @@ namespace nm { const DType* b = reinterpret_cast(b_elements); DType* x = reinterpret_cast(x_elements); - for (int i = 0; i < r; ++i) { x[i] = b[i]; } + for (int i = 0; i < r; ++i) { x[i] = b[i]; } for (int i = 0; i < r; ++i) { // forward substitution loop ip = pivot[i]; sum = x[ip]; @@ -335,18 +335,18 @@ namespace nm { for (int row = k + 1; row < M; ++row) { typename MagnitudeDType::type big; big = magnitude( matrix[M*row + k] ); // element below the temp pivot - + if ( big > akk ) { interchange = row; - akk = big; + akk = big; } - } + } if (interchange != k) { // check if rows need flipping DType temp; for (int col = 0; col < M; ++col) { - NM_SWAP(matrix[interchange*M + col], matrix[k*M + col], temp); + NM_SWAP(matrix[interchange*M + col], matrix[k*M + col], temp); } } @@ -360,7 +360,7 @@ namespace nm { DType pivot = matrix[k * (M + 1)]; matrix[k * (M + 1)] = (DType)(1); // set diagonal as 1 for in-place inversion - for (int col = 0; col < M; ++col) { + for (int col = 0; col < M; ++col) { // divide each element in the kth row with the pivot matrix[k*M + col] = matrix[k*M + col] / pivot; } @@ -369,7 +369,7 @@ namespace nm { if (kk == k) continue; DType dum = matrix[k + M*kk]; - matrix[k + M*kk] = (DType)(0); // prepare for inplace inversion + matrix[k + M*kk] = (DType)(0); // prepare for inplace inversion for (int col = 0; col < M; ++col) { matrix[M*kk + col] = matrix[M*kk + col] - matrix[M*k + col] * dum; } @@ -384,7 +384,7 @@ namespace nm { for (int row = 0; row < M; ++row) { NM_SWAP(matrix[row * M + row_index[k]], matrix[row * M + col_index[k]], - temp); + temp); } } } @@ -410,14 +410,14 @@ namespace nm { DType sum_of_squares, *p_row, *psubdiag, *p_a, scale, innerproduct; int i, k, col; - // For each column use a Householder transformation to zero all entries + // For each column use a Householder transformation to zero all entries // below the subdiagonal. - for (psubdiag = a + nrows, col = 0; col < nrows - 2; psubdiag += nrows + 1, + for (psubdiag = a + nrows, col = 0; col < nrows - 2; psubdiag += nrows + 1, col++) { // Calculate the signed square root of the sum of squares of the // elements below the diagonal. - for (p_a = psubdiag, sum_of_squares = 0.0, i = col + 1; i < nrows; + for (p_a = psubdiag, sum_of_squares = 0.0, i = col + 1; i < nrows; p_a += nrows, i++) { sum_of_squares += *p_a * *p_a; } @@ -447,7 +447,7 @@ namespace nm { *p_a -= u[k] * innerproduct; } } - + // Postmultiply QA by Q for (p_row = a, i = 0; i < nrows; p_row += nrows, i++) { for (innerproduct = 0.0, k = col + 1; k < nrows; k++) { @@ -465,7 +465,7 @@ namespace nm { } void raise_not_invertible_error() { - rb_raise(nm_eNotInvertibleError, + rb_raise(nm_eNotInvertibleError, "matrix must have non-zero determinant to be invertible (not getting this error does not mean matrix is invertible if you're dealing with floating points)"); } @@ -473,7 +473,7 @@ namespace nm { * Calculate the exact inverse for a dense matrix (A [elements]) of size 2 or 3. Places the result in B_elements. */ template - void inverse_exact_from_dense(const int M, const void* A_elements, + void inverse_exact_from_dense(const int M, const void* A_elements, const int lda, void* B_elements, const int ldb) { const DType* A = reinterpret_cast(A_elements); @@ -485,7 +485,7 @@ namespace nm { B[0] = A[lda+1] / det; B[1] = -A[1] / det; B[ldb] = -A[lda] / det; - B[ldb+1] = -A[0] / det; + B[ldb+1] = A[0] / det; } else if (M == 3) { // Calculate the exact determinant. @@ -510,7 +510,7 @@ namespace nm { } template - void inverse_exact_from_yale(const int M, const YALE_STORAGE* storage, + void inverse_exact_from_yale(const int M, const YALE_STORAGE* storage, const int lda, YALE_STORAGE* inverse, const int ldb) { // inverse is a clone of storage @@ -524,18 +524,18 @@ namespace nm { if (M == 2) { IType ndnz = ija[2] - ija[0]; - if (ndnz == 2) { - det = a[0] * a[1] - a[col_pos] * a[col_pos+1]; + if (ndnz == 2) { + det = a[0] * a[1] - a[col_pos] * a[col_pos+1]; } else { det = a[0] * a[1]; } if (det == 0) { raise_not_invertible_error(); } b[0] = a[1] / det; - b[1] = -a[0] / det; - if (ndnz == 2) { + b[1] = a[0] / det; + if (ndnz == 2) { b[col_pos] = -a[col_pos] / det; b[col_pos+1] = -a[col_pos+1] / det; } - else if (ndnz == 1) { + else if (ndnz == 1) { b[col_pos] = -a[col_pos] / det; } @@ -561,7 +561,7 @@ namespace nm { rb_raise(rb_eArgError, "some value in IJA is incorrect!"); } } - det = + det = A[0] * A[lda+1] * A[2*lda+2] + A[1] * A[lda+2] * A[2*lda] + A[2] * A[lda] * A[2*lda+1] - A[0] * A[lda+2] * A[2*lda+1] - A[1] * A[lda] * A[2*lda+2] - A[2] * A[lda+1] * A[2*lda]; if (det == 0) { raise_not_invertible_error(); } @@ -1267,9 +1267,9 @@ static VALUE nm_clapack_laswp(VALUE self, VALUE n, VALUE a, VALUE lda, VALUE k1, /* * C accessor for calculating an exact determinant. Dense matrix version. */ -void nm_math_det_exact_from_dense(const int M, const void* elements, const int lda, +void nm_math_det_exact_from_dense(const int M, const void* elements, const int lda, nm::dtype_t dtype, void* result) { - NAMED_DTYPE_TEMPLATE_TABLE(ttable, nm::math::det_exact_from_dense, void, const int M, + NAMED_DTYPE_TEMPLATE_TABLE(ttable, nm::math::det_exact_from_dense, void, const int M, const void* A_elements, const int lda, void* result_arg); ttable[dtype](M, elements, lda, result); @@ -1278,27 +1278,27 @@ void nm_math_det_exact_from_dense(const int M, const void* elements, const int l /* * C accessor for calculating an exact determinant. Yale matrix version. */ -void nm_math_det_exact_from_yale(const int M, const YALE_STORAGE* storage, const int lda, +void nm_math_det_exact_from_yale(const int M, const YALE_STORAGE* storage, const int lda, nm::dtype_t dtype, void* result) { - NAMED_DTYPE_TEMPLATE_TABLE(ttable, nm::math::det_exact_from_yale, void, const int M, + NAMED_DTYPE_TEMPLATE_TABLE(ttable, nm::math::det_exact_from_yale, void, const int M, const YALE_STORAGE* storage, const int lda, void* result_arg); ttable[dtype](M, storage, lda, result); } /* - * C accessor for solving a system of linear equations. + * C accessor for solving a system of linear equations. */ void nm_math_solve(VALUE lu, VALUE b, VALUE x, VALUE ipiv) { int* pivot = new int[RARRAY_LEN(ipiv)]; for (int i = 0; i < RARRAY_LEN(ipiv); ++i) { - pivot[i] = FIX2INT(rb_ary_entry(ipiv, i)); + pivot[i] = FIX2INT(rb_ary_entry(ipiv, i)); } NAMED_DTYPE_TEMPLATE_TABLE(ttable, nm::math::solve, void, const int, const void*, const void*, void*, const int*); - ttable[NM_DTYPE(x)](NM_SHAPE0(b), NM_STORAGE_DENSE(lu)->elements, + ttable[NM_DTYPE(x)](NM_SHAPE0(b), NM_STORAGE_DENSE(lu)->elements, NM_STORAGE_DENSE(b)->elements, NM_STORAGE_DENSE(x)->elements, pivot); } @@ -1313,7 +1313,7 @@ void nm_math_hessenberg(VALUE a) { NULL, NULL, // does not support Complex NULL // no support for Ruby Object }; - + ttable[NM_DTYPE(a)](NM_SHAPE0(a), NM_STORAGE_DENSE(a)->elements); } /* @@ -1328,10 +1328,10 @@ void nm_math_inverse(const int M, void* a_elements, nm::dtype_t dtype) { /* * C accessor for calculating an exact inverse. Dense matrix version. */ -void nm_math_inverse_exact_from_dense(const int M, const void* A_elements, +void nm_math_inverse_exact_from_dense(const int M, const void* A_elements, const int lda, void* B_elements, const int ldb, nm::dtype_t dtype) { - NAMED_DTYPE_TEMPLATE_TABLE(ttable, nm::math::inverse_exact_from_dense, void, + NAMED_DTYPE_TEMPLATE_TABLE(ttable, nm::math::inverse_exact_from_dense, void, const int, const void*, const int, void*, const int); ttable[dtype](M, A_elements, lda, B_elements, ldb); @@ -1340,10 +1340,10 @@ void nm_math_inverse_exact_from_dense(const int M, const void* A_elements, /* * C accessor for calculating an exact inverse. Yale matrix version. */ -void nm_math_inverse_exact_from_yale(const int M, const YALE_STORAGE* storage, +void nm_math_inverse_exact_from_yale(const int M, const YALE_STORAGE* storage, const int lda, YALE_STORAGE* inverse, const int ldb, nm::dtype_t dtype) { - NAMED_DTYPE_TEMPLATE_TABLE(ttable, nm::math::inverse_exact_from_yale, void, + NAMED_DTYPE_TEMPLATE_TABLE(ttable, nm::math::inverse_exact_from_yale, void, const int, const YALE_STORAGE*, const int, YALE_STORAGE*, const int); ttable[dtype](M, storage, lda, inverse, ldb); diff --git a/spec/00_nmatrix_spec.rb b/spec/00_nmatrix_spec.rb index 9ce6602b..09da20c4 100644 --- a/spec/00_nmatrix_spec.rb +++ b/spec/00_nmatrix_spec.rb @@ -73,7 +73,7 @@ e = NMatrix.new(2, [3,1,2,1], stype: :dense, dtype: :int64) inversed = e.method(:__inverse_exact__).call(e.clone, 2, 2) - f = NMatrix.new(2, [1,-1,-2,-3], stype: :dense, dtype: :int64) + f = NMatrix.new(2, [1,-1,-2,3], stype: :dense, dtype: :int64) expect(inversed).to eq(f) end @@ -91,7 +91,7 @@ e = NMatrix.new(2, [3,1,2,1], stype: :yale, dtype: :int64) inversed = e.method(:__inverse_exact__).call(e.clone, 2, 2) - f = NMatrix.new(2, [1,-1,-2,-3], stype: :yale, dtype: :int64) + f = NMatrix.new(2, [1,-1,-2,3], stype: :yale, dtype: :int64) expect(inversed).to eq(f) end @@ -104,7 +104,7 @@ c = NMatrix.new(2, [3,1,2,1], stype: :list, dtype: :int64) inversed = c.method(:__inverse_exact__).call(c.clone, 2, 2) - d = NMatrix.new(2, [1,-1,-2,-3], stype: :list, dtype: :int64) + d = NMatrix.new(2, [1,-1,-2,3], stype: :list, dtype: :int64) expect(inversed).to eq(d) end From 72d0b546b15b53a6af1220cc0da0f64d7179e835 Mon Sep 17 00:00:00 2001 From: npriyadarshi Date: Wed, 15 Mar 2017 20:48:59 +0530 Subject: [PATCH 2/4] added inverse_exact# in lib/nmatrix/math.rb --- lib/nmatrix/math.rb | 56 ++++++++++++++++++++++++++++++++++++++++++++- spec/math_spec.rb | 26 +++++++++++++++++---- 2 files changed, 76 insertions(+), 6 deletions(-) diff --git a/lib/nmatrix/math.rb b/lib/nmatrix/math.rb index d2f62710..fb91ff35 100644 --- a/lib/nmatrix/math.rb +++ b/lib/nmatrix/math.rb @@ -112,6 +112,60 @@ def invert end alias :inverse :invert + # call-seq: + # invert_exact! -> NMatrix + # + # Calulates inverse_exact of a matrix of size 2 or 3. + # Only works on dense matrices. + # + # * *Raises* : + # - +StorageTypeError+ -> only implemented on dense matrices. + # - +DataTypeError+ -> cannot invert an integer matrix in-place. + # - +NotImplementedError+ -> cannot find exact inverse of matrix with size greater than 3 # + def invert_exact! + raise(StorageTypeError, "invert only works on dense matrices currently") unless self.dense? + raise(ShapeError, "Cannot invert non-square matrix") unless self.dim == 2 && self.shape[0] == self.shape[1] + raise(DataTypeError, "Cannot invert an integer matrix in-place") if self.integer_dtype? + #No internal implementation of getri, so use this other function + n = self.shape[0] + if n>3 + raise(NotImplementedError, "Cannot find exact inverse of matrix of size greater than 3") + else + clond=self.clone + __inverse_exact__(clond, n, n) + end + end + + # + # call-seq: + # invert_exact -> NMatrix + # + # Make a copy of the matrix, then invert using exact_inverse + # + # * *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. + # + # * *Raises* : + # - +StorageTypeError+ -> only implemented on dense matrices. + # - +ShapeError+ -> matrix must be square. + # - +NotImplementedError+ -> cannot find exact inverse of matrix with size greater than 3 + # + def invert_exact + #write this in terms of invert_exact! so plugins will only have to overwrite + #invert_exact! and not invert_exact + if self.integer_dtype? + cloned = self.cast(dtype: :float64) + cloned.invert_exact! + else + cloned = self.clone + cloned.invert_exact! + end + end + alias :inverse_exact :invert_exact + + # # call-seq: @@ -1458,4 +1512,4 @@ def dtype_for_floor_or_ceil self.__dense_map__ { |l| l.send(op,rhs) } end end -end \ No newline at end of file +end diff --git a/spec/math_spec.rb b/spec/math_spec.rb index 14fd8e83..9495ec61 100644 --- a/spec/math_spec.rb +++ b/spec/math_spec.rb @@ -462,9 +462,9 @@ it "should correctly invert a matrix in place (bang)" do pending("not yet implemented for :object dtype") if dtype == :object - 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, @@ -495,6 +495,22 @@ expect(a.invert).to be_within(err).of(b) end + + it "should correctly find exact inverse" do + pending("not yet implemented for NMatrix-JRuby") if jruby? + a = NMatrix.new(:dense, 3, [1,2,3,0,1,4,5,6,0], dtype) + b = NMatrix.new(:dense, 3, [-24,18,5,20,-15,-4,-5,4,1], dtype) + + expect(a.invert_exact).to be_within(err).of(b) + end + + it "should correctly find exact inverse" do + pending("not yet implemented for NMatrix-JRuby") if jruby? + a = NMatrix.new(:dense, 2, [1,3,3,8,], dtype) + b = NMatrix.new(:dense, 2, [-8,3,3,-1], dtype) + + expect(a.invert_exact).to be_within(err).of(b) + end end end @@ -1148,7 +1164,7 @@ 3, 4, 5, 6, 7, 8], stype: stype, dtype: dtype) end - + it "scales the matrix by a given factor and return the result" do pending("not yet implemented for :object dtype") if dtype == :object if integer_dtype? dtype @@ -1160,7 +1176,7 @@ 12, 14, 16], stype: stype, dtype: dtype)) end end - + it "scales the matrix in place by a given factor" do pending("not yet implemented for :object dtype") if dtype == :object if dtype == :int8 From 4e7c5c7c402fe9d4a7564e44c5f3ce826416cd2b Mon Sep 17 00:00:00 2001 From: npriyadarshi Date: Wed, 15 Mar 2017 22:45:42 +0530 Subject: [PATCH 3/4] Removed the StorageTypeError --- lib/nmatrix/math.rb | 2 -- 1 file changed, 2 deletions(-) diff --git a/lib/nmatrix/math.rb b/lib/nmatrix/math.rb index fb91ff35..1fc55603 100644 --- a/lib/nmatrix/math.rb +++ b/lib/nmatrix/math.rb @@ -119,11 +119,9 @@ def invert # Only works on dense matrices. # # * *Raises* : - # - +StorageTypeError+ -> only implemented on dense matrices. # - +DataTypeError+ -> cannot invert an integer matrix in-place. # - +NotImplementedError+ -> cannot find exact inverse of matrix with size greater than 3 # def invert_exact! - raise(StorageTypeError, "invert only works on dense matrices currently") unless self.dense? raise(ShapeError, "Cannot invert non-square matrix") unless self.dim == 2 && self.shape[0] == self.shape[1] raise(DataTypeError, "Cannot invert an integer matrix in-place") if self.integer_dtype? #No internal implementation of getri, so use this other function From 3e1b1d7a882d74256c0ccdb387fde27629d03f35 Mon Sep 17 00:00:00 2001 From: npriyadarshi Date: Wed, 22 Mar 2017 02:16:21 +0530 Subject: [PATCH 4/4] Changed the name of the method to exact_inverse! --- lib/nmatrix/math.rb | 18 +++++++++--------- spec/math_spec.rb | 8 ++++---- 2 files changed, 13 insertions(+), 13 deletions(-) diff --git a/lib/nmatrix/math.rb b/lib/nmatrix/math.rb index 1fc55603..1883da9d 100644 --- a/lib/nmatrix/math.rb +++ b/lib/nmatrix/math.rb @@ -113,7 +113,7 @@ def invert alias :inverse :invert # call-seq: - # invert_exact! -> NMatrix + # exact_inverse! -> NMatrix # # Calulates inverse_exact of a matrix of size 2 or 3. # Only works on dense matrices. @@ -121,7 +121,7 @@ def invert # * *Raises* : # - +DataTypeError+ -> cannot invert an integer matrix in-place. # - +NotImplementedError+ -> cannot find exact inverse of matrix with size greater than 3 # - def invert_exact! + def exact_inverse! raise(ShapeError, "Cannot invert non-square matrix") unless self.dim == 2 && self.shape[0] == self.shape[1] raise(DataTypeError, "Cannot invert an integer matrix in-place") if self.integer_dtype? #No internal implementation of getri, so use this other function @@ -136,7 +136,7 @@ def invert_exact! # # call-seq: - # invert_exact -> NMatrix + # exact_inverse -> NMatrix # # Make a copy of the matrix, then invert using exact_inverse # @@ -150,18 +150,18 @@ def invert_exact! # - +ShapeError+ -> matrix must be square. # - +NotImplementedError+ -> cannot find exact inverse of matrix with size greater than 3 # - def invert_exact - #write this in terms of invert_exact! so plugins will only have to overwrite - #invert_exact! and not invert_exact + def exact_inverse + #write this in terms of exact_inverse! so plugins will only have to overwrite + #exact_inverse! and not exact_inverse if self.integer_dtype? cloned = self.cast(dtype: :float64) - cloned.invert_exact! + cloned.exact_inverse! else cloned = self.clone - cloned.invert_exact! + cloned.exact_inverse! end end - alias :inverse_exact :invert_exact + alias :invert_exactly :exact_inverse diff --git a/spec/math_spec.rb b/spec/math_spec.rb index 9495ec61..4276de6f 100644 --- a/spec/math_spec.rb +++ b/spec/math_spec.rb @@ -501,16 +501,16 @@ a = NMatrix.new(:dense, 3, [1,2,3,0,1,4,5,6,0], dtype) b = NMatrix.new(:dense, 3, [-24,18,5,20,-15,-4,-5,4,1], dtype) - expect(a.invert_exact).to be_within(err).of(b) + expect(a.exact_inverse).to be_within(err).of(b) end it "should correctly find exact inverse" do pending("not yet implemented for NMatrix-JRuby") if jruby? - a = NMatrix.new(:dense, 2, [1,3,3,8,], dtype) + a = NMatrix.new(:dense, 2, [1,3,3,8], dtype) b = NMatrix.new(:dense, 2, [-8,3,3,-1], dtype) - expect(a.invert_exact).to be_within(err).of(b) - end + expect(a.exact_inverse).to be_within(err).of(b) + end end end