From fc46b8622e7676fff8d5129b2164a9847ced418b Mon Sep 17 00:00:00 2001 From: Chris Date: Tue, 9 Jul 2019 18:16:12 -0700 Subject: [PATCH] Finished with complex arithematic --- src/apatite/linear_algebra/matrix.cr | 118 ++++++++++++++++-- .../matrix/eigenvalue_decomposition.cr | 17 +++ .../matrix/lup_decomposition.cr | 17 +++ 3 files changed, 145 insertions(+), 7 deletions(-) create mode 100644 src/apatite/linear_algebra/matrix/eigenvalue_decomposition.cr create mode 100644 src/apatite/linear_algebra/matrix/lup_decomposition.cr diff --git a/src/apatite/linear_algebra/matrix.cr b/src/apatite/linear_algebra/matrix.cr index 7ce96d8..aecda84 100644 --- a/src/apatite/linear_algebra/matrix.cr +++ b/src/apatite/linear_algebra/matrix.cr @@ -1,5 +1,6 @@ require "json" -require "./vector" +require "./matrix/eigenvalue_decomposition" +require "./matrix/lup_decomposition" module Apatite::LinearAlgebra class Matrix(T) @@ -1064,12 +1065,6 @@ module Apatite::LinearAlgebra # MATRIX FUNCTIONS -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- #++ - # Attempt to coerce the elements in the matrix to another type. - def coerce(klass) - rws = @rows.map { |r| r.map { |i| klass.new(i) } } - Matrix.rows(rws) - end - # Returns the determinant of the matrix. # # Beware that using Float values can yield erroneous results @@ -1262,6 +1257,115 @@ module Apatite::LinearAlgebra self.class.vstack(self, *matrices) end + #-- + # DECOMPOSITIONS -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= + #++ + + # Returns the Eigensystem of the matrix + # See `EigenvalueDecomposition`. + # + # NOTE: Not working yet + # + # ``` + # m = Matrix[[1, 2], [3, 4]] + # v, d, v_inv = m.eigensystem + # d.diagonal? # => true + # v.inv == v_inv # => true + # (v * d * v_inv).round(5) == m # => true + # ``` + def eigensystem + EigenvalueDecomposition.new(self) + end + + # Returns the LUP decomposition of the matrix + # See +LUPDecomposition+. + # + # NOTE: Not working yet + # + # ``` + # a = Matrix[[1, 2], [3, 4]] + # l, u, p = a.lup + # l.lower_triangular? # => true + # u.upper_triangular? # => true + # p.permutation? # => true + # l * u == p * a # => true + # a.lup.solve([2, 5]) # => Vector[(1/1), (1/2)] + # ``` + def lup + LUPDecomposition.new(self) + end + + #-- + # COMPLEX ARITHMETIC -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= + #++ + + # Returns the conjugate of the matrix. + # + # ``` + # Matrix[[Complex(1,2), Complex(0,1), 0], [1, 2, 3]] + # # => 1+2i i 0 + # # 1 2 3 + # Matrix[[Complex(1,2), Complex(0,1), 0], [1, 2, 3]].conj + # # => 1-2i -i 0 + # # 1 2 3 + # ``` + def conj + raise ArgumentError.new("Matrix#conj only works with real matrices (i.e. Matrix(Complex))") unless real? + map(&.conj) + end + + # Returns the imaginary part of the matrix. + # + # ``` + # Matrix[[Complex(1,2), Complex(0,1), 0], [1, 2, 3]] + # # => [ 1+2i, i, 0, + # # 1, 2, 3 ] + # Matrix[[Complex(1,2), Complex(0,1), 0], [1, 2, 3]].imag + # # => [ 2i, i, 0, + # # 0, 0, 0 ] + # ``` + def imag + raise ArgumentError.new("Matrix#imag only works with real matrices (i.e. Matrix(Complex))") unless real? + map(&.imag) + end + + # Returns the real part of the matrix. + # + # ``` + # Matrix[[Complex(1,2), Complex(0,1), 0], [1, 2, 3]] + # # => [ 1+2i, i, 0, + # # 1, 2, 3 ] + # Matrix[[Complex(1,2), Complex(0,1), 0], [1, 2, 3]].real + # # => [ 1, 0, 0, + # # 1, 2, 3 ] + # ``` + def real + raise ArgumentError.new("Matrix#real only works with real matrices (i.e. Matrix(Complex))") unless real? + map(&.real) + end + + # Returns an array containing matrices corresponding to the real and imaginary + # parts of the matrix + # + # ``` + # m.rect == [m.real, m.imag] + # # ==> true for all matrices m + # ``` + def rect + raise ArgumentError.new("Matrix#real only works with real matrices (i.e. Matrix(Complex))") unless real? + [real, imag] + end + + #-- + # CONVERTING -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- + #++ + + # Attempt to coerce the elements in the matrix to another type. + def coerce(klass) + rws = @rows.map { |r| r.map { |i| klass.new(i) } } + Matrix.rows(rws) + end + def to_s(io) if empty? "Matrix.empty(#{row_count}, #{column_count})" diff --git a/src/apatite/linear_algebra/matrix/eigenvalue_decomposition.cr b/src/apatite/linear_algebra/matrix/eigenvalue_decomposition.cr new file mode 100644 index 0000000..0402614 --- /dev/null +++ b/src/apatite/linear_algebra/matrix/eigenvalue_decomposition.cr @@ -0,0 +1,17 @@ +module Apatite + class Matrix(T) + # Eigenvalues and eigenvectors of a real matrix. + # + # Computes the eigenvalues and eigenvectors of a matrix A. + # + # If A is diagonalizable, this provides matrices V and D + # such that A = V*D*V.inv, where D is the diagonal matrix with entries + # equal to the eigenvalues and V is formed by the eigenvectors. + # + # If A is symmetric, then V is orthogonal and thus A = V*D*V.t + class EigenvalueDecomposition + def initialize(matrix) + end + end + end +end diff --git a/src/apatite/linear_algebra/matrix/lup_decomposition.cr b/src/apatite/linear_algebra/matrix/lup_decomposition.cr new file mode 100644 index 0000000..5404eaf --- /dev/null +++ b/src/apatite/linear_algebra/matrix/lup_decomposition.cr @@ -0,0 +1,17 @@ +module Apatite + class Matrix(T) + # For an m-by-n matrix A with m >= n, the LU decomposition is an m-by-n + # unit lower triangular matrix L, an n-by-n upper triangular matrix U, + # and a m-by-m permutation matrix P so that L*U = P*A. + # If m < n, then L is m-by-m and U is m-by-n. + # + # The LUP decomposition with pivoting always exists, even if the matrix is + # singular, so the constructor will never fail. The primary use of the + # LU decomposition is in the solution of square systems of simultaneous + # linear equations. This will fail if singular? returns true. + class LupDecomposition + def initialize(matrix) + end + end + end +end