1# encoding: utf-8 2# 3# = matrix.rb 4# 5# An implementation of Matrix and Vector classes. 6# 7# See classes Matrix and Vector for documentation. 8# 9# Current Maintainer:: Marc-André Lafortune 10# Original Author:: Keiju ISHITSUKA 11# Original Documentation:: Gavin Sinclair (sourced from <i>Ruby in a Nutshell</i> (Matsumoto, O'Reilly)) 12## 13 14require "e2mmap.rb" 15 16module ExceptionForMatrix # :nodoc: 17 extend Exception2MessageMapper 18 def_e2message(TypeError, "wrong argument type %s (expected %s)") 19 def_e2message(ArgumentError, "Wrong # of arguments(%d for %d)") 20 21 def_exception("ErrDimensionMismatch", "\#{self.name} dimension mismatch") 22 def_exception("ErrNotRegular", "Not Regular Matrix") 23 def_exception("ErrOperationNotDefined", "Operation(%s) can\\'t be defined: %s op %s") 24 def_exception("ErrOperationNotImplemented", "Sorry, Operation(%s) not implemented: %s op %s") 25end 26 27# 28# The +Matrix+ class represents a mathematical matrix. It provides methods for creating 29# matrices, operating on them arithmetically and algebraically, 30# and determining their mathematical properties (trace, rank, inverse, determinant). 31# 32# == Method Catalogue 33# 34# To create a matrix: 35# * Matrix[*rows] 36# * Matrix.[](*rows) 37# * Matrix.rows(rows, copy = true) 38# * Matrix.columns(columns) 39# * Matrix.build(row_count, column_count, &block) 40# * Matrix.diagonal(*values) 41# * Matrix.scalar(n, value) 42# * Matrix.identity(n) 43# * Matrix.unit(n) 44# * Matrix.I(n) 45# * Matrix.zero(n) 46# * Matrix.row_vector(row) 47# * Matrix.column_vector(column) 48# 49# To access Matrix elements/columns/rows/submatrices/properties: 50# * #[](i, j) 51# * #row_count (row_size) 52# * #column_count (column_size) 53# * #row(i) 54# * #column(j) 55# * #collect 56# * #map 57# * #each 58# * #each_with_index 59# * #find_index 60# * #minor(*param) 61# 62# Properties of a matrix: 63# * #diagonal? 64# * #empty? 65# * #hermitian? 66# * #lower_triangular? 67# * #normal? 68# * #orthogonal? 69# * #permutation? 70# * #real? 71# * #regular? 72# * #singular? 73# * #square? 74# * #symmetric? 75# * #unitary? 76# * #upper_triangular? 77# * #zero? 78# 79# Matrix arithmetic: 80# * #*(m) 81# * #+(m) 82# * #-(m) 83# * #/(m) 84# * #inverse 85# * #inv 86# * #** 87# 88# Matrix functions: 89# * #determinant 90# * #det 91# * #rank 92# * #round 93# * #trace 94# * #tr 95# * #transpose 96# * #t 97# 98# Matrix decompositions: 99# * #eigen 100# * #eigensystem 101# * #lup 102# * #lup_decomposition 103# 104# Complex arithmetic: 105# * conj 106# * conjugate 107# * imag 108# * imaginary 109# * real 110# * rect 111# * rectangular 112# 113# Conversion to other data types: 114# * #coerce(other) 115# * #row_vectors 116# * #column_vectors 117# * #to_a 118# 119# String representations: 120# * #to_s 121# * #inspect 122# 123class Matrix 124 include Enumerable 125 include ExceptionForMatrix 126 autoload :EigenvalueDecomposition, "matrix/eigenvalue_decomposition" 127 autoload :LUPDecomposition, "matrix/lup_decomposition" 128 129 # instance creations 130 private_class_method :new 131 attr_reader :rows 132 protected :rows 133 134 # 135 # Creates a matrix where each argument is a row. 136 # Matrix[ [25, 93], [-1, 66] ] 137 # => 25 93 138 # -1 66 139 # 140 def Matrix.[](*rows) 141 rows(rows, false) 142 end 143 144 # 145 # Creates a matrix where +rows+ is an array of arrays, each of which is a row 146 # of the matrix. If the optional argument +copy+ is false, use the given 147 # arrays as the internal structure of the matrix without copying. 148 # Matrix.rows([[25, 93], [-1, 66]]) 149 # => 25 93 150 # -1 66 151 # 152 def Matrix.rows(rows, copy = true) 153 rows = convert_to_array(rows) 154 rows.map! do |row| 155 convert_to_array(row, copy) 156 end 157 size = (rows[0] || []).size 158 rows.each do |row| 159 raise ErrDimensionMismatch, "row size differs (#{row.size} should be #{size})" unless row.size == size 160 end 161 new rows, size 162 end 163 164 # 165 # Creates a matrix using +columns+ as an array of column vectors. 166 # Matrix.columns([[25, 93], [-1, 66]]) 167 # => 25 -1 168 # 93 66 169 # 170 def Matrix.columns(columns) 171 rows(columns, false).transpose 172 end 173 174 # 175 # Creates a matrix of size +row_count+ x +column_count+. 176 # It fills the values by calling the given block, 177 # passing the current row and column. 178 # Returns an enumerator if no block is given. 179 # 180 # m = Matrix.build(2, 4) {|row, col| col - row } 181 # => Matrix[[0, 1, 2, 3], [-1, 0, 1, 2]] 182 # m = Matrix.build(3) { rand } 183 # => a 3x3 matrix with random elements 184 # 185 def Matrix.build(row_count, column_count = row_count) 186 row_count = CoercionHelper.coerce_to_int(row_count) 187 column_count = CoercionHelper.coerce_to_int(column_count) 188 raise ArgumentError if row_count < 0 || column_count < 0 189 return to_enum :build, row_count, column_count unless block_given? 190 rows = Array.new(row_count) do |i| 191 Array.new(column_count) do |j| 192 yield i, j 193 end 194 end 195 new rows, column_count 196 end 197 198 # 199 # Creates a matrix where the diagonal elements are composed of +values+. 200 # Matrix.diagonal(9, 5, -3) 201 # => 9 0 0 202 # 0 5 0 203 # 0 0 -3 204 # 205 def Matrix.diagonal(*values) 206 size = values.size 207 rows = Array.new(size) {|j| 208 row = Array.new(size, 0) 209 row[j] = values[j] 210 row 211 } 212 new rows 213 end 214 215 # 216 # Creates an +n+ by +n+ diagonal matrix where each diagonal element is 217 # +value+. 218 # Matrix.scalar(2, 5) 219 # => 5 0 220 # 0 5 221 # 222 def Matrix.scalar(n, value) 223 diagonal(*Array.new(n, value)) 224 end 225 226 # 227 # Creates an +n+ by +n+ identity matrix. 228 # Matrix.identity(2) 229 # => 1 0 230 # 0 1 231 # 232 def Matrix.identity(n) 233 scalar(n, 1) 234 end 235 class << Matrix 236 alias unit identity 237 alias I identity 238 end 239 240 # 241 # Creates a zero matrix. 242 # Matrix.zero(2) 243 # => 0 0 244 # 0 0 245 # 246 def Matrix.zero(row_count, column_count = row_count) 247 rows = Array.new(row_count){Array.new(column_count, 0)} 248 new rows, column_count 249 end 250 251 # 252 # Creates a single-row matrix where the values of that row are as given in 253 # +row+. 254 # Matrix.row_vector([4,5,6]) 255 # => 4 5 6 256 # 257 def Matrix.row_vector(row) 258 row = convert_to_array(row) 259 new [row] 260 end 261 262 # 263 # Creates a single-column matrix where the values of that column are as given 264 # in +column+. 265 # Matrix.column_vector([4,5,6]) 266 # => 4 267 # 5 268 # 6 269 # 270 def Matrix.column_vector(column) 271 column = convert_to_array(column) 272 new [column].transpose, 1 273 end 274 275 # 276 # Creates a empty matrix of +row_count+ x +column_count+. 277 # At least one of +row_count+ or +column_count+ must be 0. 278 # 279 # m = Matrix.empty(2, 0) 280 # m == Matrix[ [], [] ] 281 # => true 282 # n = Matrix.empty(0, 3) 283 # n == Matrix.columns([ [], [], [] ]) 284 # => true 285 # m * n 286 # => Matrix[[0, 0, 0], [0, 0, 0]] 287 # 288 def Matrix.empty(row_count = 0, column_count = 0) 289 raise ArgumentError, "One size must be 0" if column_count != 0 && row_count != 0 290 raise ArgumentError, "Negative size" if column_count < 0 || row_count < 0 291 292 new([[]]*row_count, column_count) 293 end 294 295 # 296 # Matrix.new is private; use Matrix.rows, columns, [], etc... to create. 297 # 298 def initialize(rows, column_count = rows[0].size) 299 # No checking is done at this point. rows must be an Array of Arrays. 300 # column_count must be the size of the first row, if there is one, 301 # otherwise it *must* be specified and can be any integer >= 0 302 @rows = rows 303 @column_count = column_count 304 end 305 306 def new_matrix(rows, column_count = rows[0].size) # :nodoc: 307 self.class.send(:new, rows, column_count) # bypass privacy of Matrix.new 308 end 309 private :new_matrix 310 311 # 312 # Returns element (+i+,+j+) of the matrix. That is: row +i+, column +j+. 313 # 314 def [](i, j) 315 @rows.fetch(i){return nil}[j] 316 end 317 alias element [] 318 alias component [] 319 320 def []=(i, j, v) 321 @rows[i][j] = v 322 end 323 alias set_element []= 324 alias set_component []= 325 private :[]=, :set_element, :set_component 326 327 # 328 # Returns the number of rows. 329 # 330 def row_count 331 @rows.size 332 end 333 334 alias_method :row_size, :row_count 335 # 336 # Returns the number of columns. 337 # 338 attr_reader :column_count 339 alias_method :column_size, :column_count 340 341 # 342 # Returns row vector number +i+ of the matrix as a Vector (starting at 0 like 343 # an array). When a block is given, the elements of that vector are iterated. 344 # 345 def row(i, &block) # :yield: e 346 if block_given? 347 @rows.fetch(i){return self}.each(&block) 348 self 349 else 350 Vector.elements(@rows.fetch(i){return nil}) 351 end 352 end 353 354 # 355 # Returns column vector number +j+ of the matrix as a Vector (starting at 0 356 # like an array). When a block is given, the elements of that vector are 357 # iterated. 358 # 359 def column(j) # :yield: e 360 if block_given? 361 return self if j >= column_count || j < -column_count 362 row_count.times do |i| 363 yield @rows[i][j] 364 end 365 self 366 else 367 return nil if j >= column_count || j < -column_count 368 col = Array.new(row_count) {|i| 369 @rows[i][j] 370 } 371 Vector.elements(col, false) 372 end 373 end 374 375 # 376 # Returns a matrix that is the result of iteration of the given block over all 377 # elements of the matrix. 378 # Matrix[ [1,2], [3,4] ].collect { |e| e**2 } 379 # => 1 4 380 # 9 16 381 # 382 def collect(&block) # :yield: e 383 return to_enum(:collect) unless block_given? 384 rows = @rows.collect{|row| row.collect(&block)} 385 new_matrix rows, column_count 386 end 387 alias map collect 388 389 # 390 # Yields all elements of the matrix, starting with those of the first row, 391 # or returns an Enumerator is no block given. 392 # Elements can be restricted by passing an argument: 393 # * :all (default): yields all elements 394 # * :diagonal: yields only elements on the diagonal 395 # * :off_diagonal: yields all elements except on the diagonal 396 # * :lower: yields only elements on or below the diagonal 397 # * :strict_lower: yields only elements below the diagonal 398 # * :strict_upper: yields only elements above the diagonal 399 # * :upper: yields only elements on or above the diagonal 400 # 401 # Matrix[ [1,2], [3,4] ].each { |e| puts e } 402 # # => prints the numbers 1 to 4 403 # Matrix[ [1,2], [3,4] ].each(:strict_lower).to_a # => [3] 404 # 405 def each(which = :all) # :yield: e 406 return to_enum :each, which unless block_given? 407 last = column_count - 1 408 case which 409 when :all 410 block = Proc.new 411 @rows.each do |row| 412 row.each(&block) 413 end 414 when :diagonal 415 @rows.each_with_index do |row, row_index| 416 yield row.fetch(row_index){return self} 417 end 418 when :off_diagonal 419 @rows.each_with_index do |row, row_index| 420 column_count.times do |col_index| 421 yield row[col_index] unless row_index == col_index 422 end 423 end 424 when :lower 425 @rows.each_with_index do |row, row_index| 426 0.upto([row_index, last].min) do |col_index| 427 yield row[col_index] 428 end 429 end 430 when :strict_lower 431 @rows.each_with_index do |row, row_index| 432 [row_index, column_count].min.times do |col_index| 433 yield row[col_index] 434 end 435 end 436 when :strict_upper 437 @rows.each_with_index do |row, row_index| 438 (row_index+1).upto(last) do |col_index| 439 yield row[col_index] 440 end 441 end 442 when :upper 443 @rows.each_with_index do |row, row_index| 444 row_index.upto(last) do |col_index| 445 yield row[col_index] 446 end 447 end 448 else 449 raise ArgumentError, "expected #{which.inspect} to be one of :all, :diagonal, :off_diagonal, :lower, :strict_lower, :strict_upper or :upper" 450 end 451 self 452 end 453 454 # 455 # Same as #each, but the row index and column index in addition to the element 456 # 457 # Matrix[ [1,2], [3,4] ].each_with_index do |e, row, col| 458 # puts "#{e} at #{row}, #{col}" 459 # end 460 # # => Prints: 461 # # 1 at 0, 0 462 # # 2 at 0, 1 463 # # 3 at 1, 0 464 # # 4 at 1, 1 465 # 466 def each_with_index(which = :all) # :yield: e, row, column 467 return to_enum :each_with_index, which unless block_given? 468 last = column_count - 1 469 case which 470 when :all 471 @rows.each_with_index do |row, row_index| 472 row.each_with_index do |e, col_index| 473 yield e, row_index, col_index 474 end 475 end 476 when :diagonal 477 @rows.each_with_index do |row, row_index| 478 yield row.fetch(row_index){return self}, row_index, row_index 479 end 480 when :off_diagonal 481 @rows.each_with_index do |row, row_index| 482 column_count.times do |col_index| 483 yield row[col_index], row_index, col_index unless row_index == col_index 484 end 485 end 486 when :lower 487 @rows.each_with_index do |row, row_index| 488 0.upto([row_index, last].min) do |col_index| 489 yield row[col_index], row_index, col_index 490 end 491 end 492 when :strict_lower 493 @rows.each_with_index do |row, row_index| 494 [row_index, column_count].min.times do |col_index| 495 yield row[col_index], row_index, col_index 496 end 497 end 498 when :strict_upper 499 @rows.each_with_index do |row, row_index| 500 (row_index+1).upto(last) do |col_index| 501 yield row[col_index], row_index, col_index 502 end 503 end 504 when :upper 505 @rows.each_with_index do |row, row_index| 506 row_index.upto(last) do |col_index| 507 yield row[col_index], row_index, col_index 508 end 509 end 510 else 511 raise ArgumentError, "expected #{which.inspect} to be one of :all, :diagonal, :off_diagonal, :lower, :strict_lower, :strict_upper or :upper" 512 end 513 self 514 end 515 516 SELECTORS = {all: true, diagonal: true, off_diagonal: true, lower: true, strict_lower: true, strict_upper: true, upper: true}.freeze 517 # 518 # :call-seq: 519 # index(value, selector = :all) -> [row, column] 520 # index(selector = :all){ block } -> [row, column] 521 # index(selector = :all) -> an_enumerator 522 # 523 # The index method is specialized to return the index as [row, column] 524 # It also accepts an optional +selector+ argument, see #each for details. 525 # 526 # Matrix[ [1,2], [3,4] ].index(&:even?) # => [0, 1] 527 # Matrix[ [1,1], [1,1] ].index(1, :strict_lower) # => [1, 0] 528 # 529 def index(*args) 530 raise ArgumentError, "wrong number of arguments(#{args.size} for 0-2)" if args.size > 2 531 which = (args.size == 2 || SELECTORS.include?(args.last)) ? args.pop : :all 532 return to_enum :find_index, which, *args unless block_given? || args.size == 1 533 if args.size == 1 534 value = args.first 535 each_with_index(which) do |e, row_index, col_index| 536 return row_index, col_index if e == value 537 end 538 else 539 each_with_index(which) do |e, row_index, col_index| 540 return row_index, col_index if yield e 541 end 542 end 543 nil 544 end 545 alias_method :find_index, :index 546 # 547 # Returns a section of the matrix. The parameters are either: 548 # * start_row, nrows, start_col, ncols; OR 549 # * row_range, col_range 550 # 551 # Matrix.diagonal(9, 5, -3).minor(0..1, 0..2) 552 # => 9 0 0 553 # 0 5 0 554 # 555 # Like Array#[], negative indices count backward from the end of the 556 # row or column (-1 is the last element). Returns nil if the starting 557 # row or column is greater than row_count or column_count respectively. 558 # 559 def minor(*param) 560 case param.size 561 when 2 562 row_range, col_range = param 563 from_row = row_range.first 564 from_row += row_count if from_row < 0 565 to_row = row_range.end 566 to_row += row_count if to_row < 0 567 to_row += 1 unless row_range.exclude_end? 568 size_row = to_row - from_row 569 570 from_col = col_range.first 571 from_col += column_count if from_col < 0 572 to_col = col_range.end 573 to_col += column_count if to_col < 0 574 to_col += 1 unless col_range.exclude_end? 575 size_col = to_col - from_col 576 when 4 577 from_row, size_row, from_col, size_col = param 578 return nil if size_row < 0 || size_col < 0 579 from_row += row_count if from_row < 0 580 from_col += column_count if from_col < 0 581 else 582 raise ArgumentError, param.inspect 583 end 584 585 return nil if from_row > row_count || from_col > column_count || from_row < 0 || from_col < 0 586 rows = @rows[from_row, size_row].collect{|row| 587 row[from_col, size_col] 588 } 589 new_matrix rows, [column_count - from_col, size_col].min 590 end 591 592 #-- 593 # TESTING -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- 594 #++ 595 596 # 597 # Returns +true+ is this is a diagonal matrix. 598 # Raises an error if matrix is not square. 599 # 600 def diagonal? 601 Matrix.Raise ErrDimensionMismatch unless square? 602 each(:off_diagonal).all?(&:zero?) 603 end 604 605 # 606 # Returns +true+ if this is an empty matrix, i.e. if the number of rows 607 # or the number of columns is 0. 608 # 609 def empty? 610 column_count == 0 || row_count == 0 611 end 612 613 # 614 # Returns +true+ is this is an hermitian matrix. 615 # Raises an error if matrix is not square. 616 # 617 def hermitian? 618 Matrix.Raise ErrDimensionMismatch unless square? 619 each_with_index(:upper).all? do |e, row, col| 620 e == rows[col][row].conj 621 end 622 end 623 624 # 625 # Returns +true+ is this is a lower triangular matrix. 626 # 627 def lower_triangular? 628 each(:strict_upper).all?(&:zero?) 629 end 630 631 # 632 # Returns +true+ is this is a normal matrix. 633 # Raises an error if matrix is not square. 634 # 635 def normal? 636 Matrix.Raise ErrDimensionMismatch unless square? 637 rows.each_with_index do |row_i, i| 638 rows.each_with_index do |row_j, j| 639 s = 0 640 rows.each_with_index do |row_k, k| 641 s += row_i[k] * row_j[k].conj - row_k[i].conj * row_k[j] 642 end 643 return false unless s == 0 644 end 645 end 646 true 647 end 648 649 # 650 # Returns +true+ is this is an orthogonal matrix 651 # Raises an error if matrix is not square. 652 # 653 def orthogonal? 654 Matrix.Raise ErrDimensionMismatch unless square? 655 rows.each_with_index do |row, i| 656 column_count.times do |j| 657 s = 0 658 row_count.times do |k| 659 s += row[k] * rows[k][j] 660 end 661 return false unless s == (i == j ? 1 : 0) 662 end 663 end 664 true 665 end 666 667 # 668 # Returns +true+ is this is a permutation matrix 669 # Raises an error if matrix is not square. 670 # 671 def permutation? 672 Matrix.Raise ErrDimensionMismatch unless square? 673 cols = Array.new(column_count) 674 rows.each_with_index do |row, i| 675 found = false 676 row.each_with_index do |e, j| 677 if e == 1 678 return false if found || cols[j] 679 found = cols[j] = true 680 elsif e != 0 681 return false 682 end 683 end 684 return false unless found 685 end 686 true 687 end 688 689 # 690 # Returns +true+ if all entries of the matrix are real. 691 # 692 def real? 693 all?(&:real?) 694 end 695 696 # 697 # Returns +true+ if this is a regular (i.e. non-singular) matrix. 698 # 699 def regular? 700 not singular? 701 end 702 703 # 704 # Returns +true+ is this is a singular matrix. 705 # 706 def singular? 707 determinant == 0 708 end 709 710 # 711 # Returns +true+ is this is a square matrix. 712 # 713 def square? 714 column_count == row_count 715 end 716 717 # 718 # Returns +true+ is this is a symmetric matrix. 719 # Raises an error if matrix is not square. 720 # 721 def symmetric? 722 Matrix.Raise ErrDimensionMismatch unless square? 723 each_with_index(:strict_upper) do |e, row, col| 724 return false if e != rows[col][row] 725 end 726 true 727 end 728 729 # 730 # Returns +true+ is this is a unitary matrix 731 # Raises an error if matrix is not square. 732 # 733 def unitary? 734 Matrix.Raise ErrDimensionMismatch unless square? 735 rows.each_with_index do |row, i| 736 column_count.times do |j| 737 s = 0 738 row_count.times do |k| 739 s += row[k].conj * rows[k][j] 740 end 741 return false unless s == (i == j ? 1 : 0) 742 end 743 end 744 true 745 end 746 747 # 748 # Returns +true+ is this is an upper triangular matrix. 749 # 750 def upper_triangular? 751 each(:strict_lower).all?(&:zero?) 752 end 753 754 # 755 # Returns +true+ is this is a matrix with only zero elements 756 # 757 def zero? 758 all?(&:zero?) 759 end 760 761 #-- 762 # OBJECT METHODS -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- 763 #++ 764 765 # 766 # Returns +true+ if and only if the two matrices contain equal elements. 767 # 768 def ==(other) 769 return false unless Matrix === other && 770 column_count == other.column_count # necessary for empty matrices 771 rows == other.rows 772 end 773 774 def eql?(other) 775 return false unless Matrix === other && 776 column_count == other.column_count # necessary for empty matrices 777 rows.eql? other.rows 778 end 779 780 # 781 # Returns a clone of the matrix, so that the contents of each do not reference 782 # identical objects. 783 # There should be no good reason to do this since Matrices are immutable. 784 # 785 def clone 786 new_matrix @rows.map(&:dup), column_count 787 end 788 789 # 790 # Returns a hash-code for the matrix. 791 # 792 def hash 793 @rows.hash 794 end 795 796 #-- 797 # ARITHMETIC -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- 798 #++ 799 800 # 801 # Matrix multiplication. 802 # Matrix[[2,4], [6,8]] * Matrix.identity(2) 803 # => 2 4 804 # 6 8 805 # 806 def *(m) # m is matrix or vector or number 807 case(m) 808 when Numeric 809 rows = @rows.collect {|row| 810 row.collect {|e| e * m } 811 } 812 return new_matrix rows, column_count 813 when Vector 814 m = self.class.column_vector(m) 815 r = self * m 816 return r.column(0) 817 when Matrix 818 Matrix.Raise ErrDimensionMismatch if column_count != m.row_count 819 820 rows = Array.new(row_count) {|i| 821 Array.new(m.column_count) {|j| 822 (0 ... column_count).inject(0) do |vij, k| 823 vij + self[i, k] * m[k, j] 824 end 825 } 826 } 827 return new_matrix rows, m.column_count 828 else 829 return apply_through_coercion(m, __method__) 830 end 831 end 832 833 # 834 # Matrix addition. 835 # Matrix.scalar(2,5) + Matrix[[1,0], [-4,7]] 836 # => 6 0 837 # -4 12 838 # 839 def +(m) 840 case m 841 when Numeric 842 Matrix.Raise ErrOperationNotDefined, "+", self.class, m.class 843 when Vector 844 m = self.class.column_vector(m) 845 when Matrix 846 else 847 return apply_through_coercion(m, __method__) 848 end 849 850 Matrix.Raise ErrDimensionMismatch unless row_count == m.row_count and column_count == m.column_count 851 852 rows = Array.new(row_count) {|i| 853 Array.new(column_count) {|j| 854 self[i, j] + m[i, j] 855 } 856 } 857 new_matrix rows, column_count 858 end 859 860 # 861 # Matrix subtraction. 862 # Matrix[[1,5], [4,2]] - Matrix[[9,3], [-4,1]] 863 # => -8 2 864 # 8 1 865 # 866 def -(m) 867 case m 868 when Numeric 869 Matrix.Raise ErrOperationNotDefined, "-", self.class, m.class 870 when Vector 871 m = self.class.column_vector(m) 872 when Matrix 873 else 874 return apply_through_coercion(m, __method__) 875 end 876 877 Matrix.Raise ErrDimensionMismatch unless row_count == m.row_count and column_count == m.column_count 878 879 rows = Array.new(row_count) {|i| 880 Array.new(column_count) {|j| 881 self[i, j] - m[i, j] 882 } 883 } 884 new_matrix rows, column_count 885 end 886 887 # 888 # Matrix division (multiplication by the inverse). 889 # Matrix[[7,6], [3,9]] / Matrix[[2,9], [3,1]] 890 # => -7 1 891 # -3 -6 892 # 893 def /(other) 894 case other 895 when Numeric 896 rows = @rows.collect {|row| 897 row.collect {|e| e / other } 898 } 899 return new_matrix rows, column_count 900 when Matrix 901 return self * other.inverse 902 else 903 return apply_through_coercion(other, __method__) 904 end 905 end 906 907 # 908 # Returns the inverse of the matrix. 909 # Matrix[[-1, -1], [0, -1]].inverse 910 # => -1 1 911 # 0 -1 912 # 913 def inverse 914 Matrix.Raise ErrDimensionMismatch unless square? 915 self.class.I(row_count).send(:inverse_from, self) 916 end 917 alias inv inverse 918 919 def inverse_from(src) # :nodoc: 920 last = row_count - 1 921 a = src.to_a 922 923 0.upto(last) do |k| 924 i = k 925 akk = a[k][k].abs 926 (k+1).upto(last) do |j| 927 v = a[j][k].abs 928 if v > akk 929 i = j 930 akk = v 931 end 932 end 933 Matrix.Raise ErrNotRegular if akk == 0 934 if i != k 935 a[i], a[k] = a[k], a[i] 936 @rows[i], @rows[k] = @rows[k], @rows[i] 937 end 938 akk = a[k][k] 939 940 0.upto(last) do |ii| 941 next if ii == k 942 q = a[ii][k].quo(akk) 943 a[ii][k] = 0 944 945 (k + 1).upto(last) do |j| 946 a[ii][j] -= a[k][j] * q 947 end 948 0.upto(last) do |j| 949 @rows[ii][j] -= @rows[k][j] * q 950 end 951 end 952 953 (k+1).upto(last) do |j| 954 a[k][j] = a[k][j].quo(akk) 955 end 956 0.upto(last) do |j| 957 @rows[k][j] = @rows[k][j].quo(akk) 958 end 959 end 960 self 961 end 962 private :inverse_from 963 964 # 965 # Matrix exponentiation. 966 # Equivalent to multiplying the matrix by itself N times. 967 # Non integer exponents will be handled by diagonalizing the matrix. 968 # 969 # Matrix[[7,6], [3,9]] ** 2 970 # => 67 96 971 # 48 99 972 # 973 def ** (other) 974 case other 975 when Integer 976 x = self 977 if other <= 0 978 x = self.inverse 979 return self.class.identity(self.column_count) if other == 0 980 other = -other 981 end 982 z = nil 983 loop do 984 z = z ? z * x : x if other[0] == 1 985 return z if (other >>= 1).zero? 986 x *= x 987 end 988 when Numeric 989 v, d, v_inv = eigensystem 990 v * self.class.diagonal(*d.each(:diagonal).map{|e| e ** other}) * v_inv 991 else 992 Matrix.Raise ErrOperationNotDefined, "**", self.class, other.class 993 end 994 end 995 996 #-- 997 # MATRIX FUNCTIONS -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- 998 #++ 999 1000 # 1001 # Returns the determinant of the matrix. 1002 # 1003 # Beware that using Float values can yield erroneous results 1004 # because of their lack of precision. 1005 # Consider using exact types like Rational or BigDecimal instead. 1006 # 1007 # Matrix[[7,6], [3,9]].determinant 1008 # => 45 1009 # 1010 def determinant 1011 Matrix.Raise ErrDimensionMismatch unless square? 1012 m = @rows 1013 case row_count 1014 # Up to 4x4, give result using Laplacian expansion by minors. 1015 # This will typically be faster, as well as giving good results 1016 # in case of Floats 1017 when 0 1018 +1 1019 when 1 1020 + m[0][0] 1021 when 2 1022 + m[0][0] * m[1][1] - m[0][1] * m[1][0] 1023 when 3 1024 m0, m1, m2 = m 1025 + m0[0] * m1[1] * m2[2] - m0[0] * m1[2] * m2[1] \ 1026 - m0[1] * m1[0] * m2[2] + m0[1] * m1[2] * m2[0] \ 1027 + m0[2] * m1[0] * m2[1] - m0[2] * m1[1] * m2[0] 1028 when 4 1029 m0, m1, m2, m3 = m 1030 + m0[0] * m1[1] * m2[2] * m3[3] - m0[0] * m1[1] * m2[3] * m3[2] \ 1031 - m0[0] * m1[2] * m2[1] * m3[3] + m0[0] * m1[2] * m2[3] * m3[1] \ 1032 + m0[0] * m1[3] * m2[1] * m3[2] - m0[0] * m1[3] * m2[2] * m3[1] \ 1033 - m0[1] * m1[0] * m2[2] * m3[3] + m0[1] * m1[0] * m2[3] * m3[2] \ 1034 + m0[1] * m1[2] * m2[0] * m3[3] - m0[1] * m1[2] * m2[3] * m3[0] \ 1035 - m0[1] * m1[3] * m2[0] * m3[2] + m0[1] * m1[3] * m2[2] * m3[0] \ 1036 + m0[2] * m1[0] * m2[1] * m3[3] - m0[2] * m1[0] * m2[3] * m3[1] \ 1037 - m0[2] * m1[1] * m2[0] * m3[3] + m0[2] * m1[1] * m2[3] * m3[0] \ 1038 + m0[2] * m1[3] * m2[0] * m3[1] - m0[2] * m1[3] * m2[1] * m3[0] \ 1039 - m0[3] * m1[0] * m2[1] * m3[2] + m0[3] * m1[0] * m2[2] * m3[1] \ 1040 + m0[3] * m1[1] * m2[0] * m3[2] - m0[3] * m1[1] * m2[2] * m3[0] \ 1041 - m0[3] * m1[2] * m2[0] * m3[1] + m0[3] * m1[2] * m2[1] * m3[0] 1042 else 1043 # For bigger matrices, use an efficient and general algorithm. 1044 # Currently, we use the Gauss-Bareiss algorithm 1045 determinant_bareiss 1046 end 1047 end 1048 alias_method :det, :determinant 1049 1050 # 1051 # Private. Use Matrix#determinant 1052 # 1053 # Returns the determinant of the matrix, using 1054 # Bareiss' multistep integer-preserving gaussian elimination. 1055 # It has the same computational cost order O(n^3) as standard Gaussian elimination. 1056 # Intermediate results are fraction free and of lower complexity. 1057 # A matrix of Integers will have thus intermediate results that are also Integers, 1058 # with smaller bignums (if any), while a matrix of Float will usually have 1059 # intermediate results with better precision. 1060 # 1061 def determinant_bareiss 1062 size = row_count 1063 last = size - 1 1064 a = to_a 1065 no_pivot = Proc.new{ return 0 } 1066 sign = +1 1067 pivot = 1 1068 size.times do |k| 1069 previous_pivot = pivot 1070 if (pivot = a[k][k]) == 0 1071 switch = (k+1 ... size).find(no_pivot) {|row| 1072 a[row][k] != 0 1073 } 1074 a[switch], a[k] = a[k], a[switch] 1075 pivot = a[k][k] 1076 sign = -sign 1077 end 1078 (k+1).upto(last) do |i| 1079 ai = a[i] 1080 (k+1).upto(last) do |j| 1081 ai[j] = (pivot * ai[j] - ai[k] * a[k][j]) / previous_pivot 1082 end 1083 end 1084 end 1085 sign * pivot 1086 end 1087 private :determinant_bareiss 1088 1089 # 1090 # deprecated; use Matrix#determinant 1091 # 1092 def determinant_e 1093 warn "#{caller(1)[0]}: warning: Matrix#determinant_e is deprecated; use #determinant" 1094 determinant 1095 end 1096 alias det_e determinant_e 1097 1098 # 1099 # Returns the rank of the matrix. 1100 # Beware that using Float values can yield erroneous results 1101 # because of their lack of precision. 1102 # Consider using exact types like Rational or BigDecimal instead. 1103 # 1104 # Matrix[[7,6], [3,9]].rank 1105 # => 2 1106 # 1107 def rank 1108 # We currently use Bareiss' multistep integer-preserving gaussian elimination 1109 # (see comments on determinant) 1110 a = to_a 1111 last_column = column_count - 1 1112 last_row = row_count - 1 1113 pivot_row = 0 1114 previous_pivot = 1 1115 0.upto(last_column) do |k| 1116 switch_row = (pivot_row .. last_row).find {|row| 1117 a[row][k] != 0 1118 } 1119 if switch_row 1120 a[switch_row], a[pivot_row] = a[pivot_row], a[switch_row] unless pivot_row == switch_row 1121 pivot = a[pivot_row][k] 1122 (pivot_row+1).upto(last_row) do |i| 1123 ai = a[i] 1124 (k+1).upto(last_column) do |j| 1125 ai[j] = (pivot * ai[j] - ai[k] * a[pivot_row][j]) / previous_pivot 1126 end 1127 end 1128 pivot_row += 1 1129 previous_pivot = pivot 1130 end 1131 end 1132 pivot_row 1133 end 1134 1135 # 1136 # deprecated; use Matrix#rank 1137 # 1138 def rank_e 1139 warn "#{caller(1)[0]}: warning: Matrix#rank_e is deprecated; use #rank" 1140 rank 1141 end 1142 1143 # Returns a matrix with entries rounded to the given precision 1144 # (see Float#round) 1145 # 1146 def round(ndigits=0) 1147 map{|e| e.round(ndigits)} 1148 end 1149 1150 # 1151 # Returns the trace (sum of diagonal elements) of the matrix. 1152 # Matrix[[7,6], [3,9]].trace 1153 # => 16 1154 # 1155 def trace 1156 Matrix.Raise ErrDimensionMismatch unless square? 1157 (0...column_count).inject(0) do |tr, i| 1158 tr + @rows[i][i] 1159 end 1160 end 1161 alias tr trace 1162 1163 # 1164 # Returns the transpose of the matrix. 1165 # Matrix[[1,2], [3,4], [5,6]] 1166 # => 1 2 1167 # 3 4 1168 # 5 6 1169 # Matrix[[1,2], [3,4], [5,6]].transpose 1170 # => 1 3 5 1171 # 2 4 6 1172 # 1173 def transpose 1174 return self.class.empty(column_count, 0) if row_count.zero? 1175 new_matrix @rows.transpose, row_count 1176 end 1177 alias t transpose 1178 1179 #-- 1180 # DECOMPOSITIONS -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= 1181 #++ 1182 1183 # 1184 # Returns the Eigensystem of the matrix; see +EigenvalueDecomposition+. 1185 # m = Matrix[[1, 2], [3, 4]] 1186 # v, d, v_inv = m.eigensystem 1187 # d.diagonal? # => true 1188 # v.inv == v_inv # => true 1189 # (v * d * v_inv).round(5) == m # => true 1190 # 1191 def eigensystem 1192 EigenvalueDecomposition.new(self) 1193 end 1194 alias eigen eigensystem 1195 1196 # 1197 # Returns the LUP decomposition of the matrix; see +LUPDecomposition+. 1198 # a = Matrix[[1, 2], [3, 4]] 1199 # l, u, p = a.lup 1200 # l.lower_triangular? # => true 1201 # u.upper_triangular? # => true 1202 # p.permutation? # => true 1203 # l * u == p * a # => true 1204 # a.lup.solve([2, 5]) # => Vector[(1/1), (1/2)] 1205 # 1206 def lup 1207 LUPDecomposition.new(self) 1208 end 1209 alias lup_decomposition lup 1210 1211 #-- 1212 # COMPLEX ARITHMETIC -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= 1213 #++ 1214 1215 # 1216 # Returns the conjugate of the matrix. 1217 # Matrix[[Complex(1,2), Complex(0,1), 0], [1, 2, 3]] 1218 # => 1+2i i 0 1219 # 1 2 3 1220 # Matrix[[Complex(1,2), Complex(0,1), 0], [1, 2, 3]].conjugate 1221 # => 1-2i -i 0 1222 # 1 2 3 1223 # 1224 def conjugate 1225 collect(&:conjugate) 1226 end 1227 alias conj conjugate 1228 1229 # 1230 # Returns the imaginary part of the matrix. 1231 # Matrix[[Complex(1,2), Complex(0,1), 0], [1, 2, 3]] 1232 # => 1+2i i 0 1233 # 1 2 3 1234 # Matrix[[Complex(1,2), Complex(0,1), 0], [1, 2, 3]].imaginary 1235 # => 2i i 0 1236 # 0 0 0 1237 # 1238 def imaginary 1239 collect(&:imaginary) 1240 end 1241 alias imag imaginary 1242 1243 # 1244 # Returns the real part of the matrix. 1245 # Matrix[[Complex(1,2), Complex(0,1), 0], [1, 2, 3]] 1246 # => 1+2i i 0 1247 # 1 2 3 1248 # Matrix[[Complex(1,2), Complex(0,1), 0], [1, 2, 3]].real 1249 # => 1 0 0 1250 # 1 2 3 1251 # 1252 def real 1253 collect(&:real) 1254 end 1255 1256 # 1257 # Returns an array containing matrices corresponding to the real and imaginary 1258 # parts of the matrix 1259 # 1260 # m.rect == [m.real, m.imag] # ==> true for all matrices m 1261 # 1262 def rect 1263 [real, imag] 1264 end 1265 alias rectangular rect 1266 1267 #-- 1268 # CONVERTING -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- 1269 #++ 1270 1271 # 1272 # The coerce method provides support for Ruby type coercion. 1273 # This coercion mechanism is used by Ruby to handle mixed-type 1274 # numeric operations: it is intended to find a compatible common 1275 # type between the two operands of the operator. 1276 # See also Numeric#coerce. 1277 # 1278 def coerce(other) 1279 case other 1280 when Numeric 1281 return Scalar.new(other), self 1282 else 1283 raise TypeError, "#{self.class} can't be coerced into #{other.class}" 1284 end 1285 end 1286 1287 # 1288 # Returns an array of the row vectors of the matrix. See Vector. 1289 # 1290 def row_vectors 1291 Array.new(row_count) {|i| 1292 row(i) 1293 } 1294 end 1295 1296 # 1297 # Returns an array of the column vectors of the matrix. See Vector. 1298 # 1299 def column_vectors 1300 Array.new(column_count) {|i| 1301 column(i) 1302 } 1303 end 1304 1305 # 1306 # Returns an array of arrays that describe the rows of the matrix. 1307 # 1308 def to_a 1309 @rows.collect(&:dup) 1310 end 1311 1312 def elements_to_f 1313 warn "#{caller(1)[0]}: warning: Matrix#elements_to_f is deprecated, use map(&:to_f)" 1314 map(&:to_f) 1315 end 1316 1317 def elements_to_i 1318 warn "#{caller(1)[0]}: warning: Matrix#elements_to_i is deprecated, use map(&:to_i)" 1319 map(&:to_i) 1320 end 1321 1322 def elements_to_r 1323 warn "#{caller(1)[0]}: warning: Matrix#elements_to_r is deprecated, use map(&:to_r)" 1324 map(&:to_r) 1325 end 1326 1327 #-- 1328 # PRINTING -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- 1329 #++ 1330 1331 # 1332 # Overrides Object#to_s 1333 # 1334 def to_s 1335 if empty? 1336 "#{self.class}.empty(#{row_count}, #{column_count})" 1337 else 1338 "#{self.class}[" + @rows.collect{|row| 1339 "[" + row.collect{|e| e.to_s}.join(", ") + "]" 1340 }.join(", ")+"]" 1341 end 1342 end 1343 1344 # 1345 # Overrides Object#inspect 1346 # 1347 def inspect 1348 if empty? 1349 "#{self.class}.empty(#{row_count}, #{column_count})" 1350 else 1351 "#{self.class}#{@rows.inspect}" 1352 end 1353 end 1354 1355 # Private helper modules 1356 1357 module ConversionHelper # :nodoc: 1358 # 1359 # Converts the obj to an Array. If copy is set to true 1360 # a copy of obj will be made if necessary. 1361 # 1362 def convert_to_array(obj, copy = false) # :nodoc: 1363 case obj 1364 when Array 1365 copy ? obj.dup : obj 1366 when Vector 1367 obj.to_a 1368 else 1369 begin 1370 converted = obj.to_ary 1371 rescue Exception => e 1372 raise TypeError, "can't convert #{obj.class} into an Array (#{e.message})" 1373 end 1374 raise TypeError, "#{obj.class}#to_ary should return an Array" unless converted.is_a? Array 1375 converted 1376 end 1377 end 1378 private :convert_to_array 1379 end 1380 1381 extend ConversionHelper 1382 1383 module CoercionHelper # :nodoc: 1384 # 1385 # Applies the operator +oper+ with argument +obj+ 1386 # through coercion of +obj+ 1387 # 1388 def apply_through_coercion(obj, oper) 1389 coercion = obj.coerce(self) 1390 raise TypeError unless coercion.is_a?(Array) && coercion.length == 2 1391 coercion[0].public_send(oper, coercion[1]) 1392 rescue 1393 raise TypeError, "#{obj.inspect} can't be coerced into #{self.class}" 1394 end 1395 private :apply_through_coercion 1396 1397 # 1398 # Helper method to coerce a value into a specific class. 1399 # Raises a TypeError if the coercion fails or the returned value 1400 # is not of the right class. 1401 # (from Rubinius) 1402 # 1403 def self.coerce_to(obj, cls, meth) # :nodoc: 1404 return obj if obj.kind_of?(cls) 1405 1406 begin 1407 ret = obj.__send__(meth) 1408 rescue Exception => e 1409 raise TypeError, "Coercion error: #{obj.inspect}.#{meth} => #{cls} failed:\n" \ 1410 "(#{e.message})" 1411 end 1412 raise TypeError, "Coercion error: obj.#{meth} did NOT return a #{cls} (was #{ret.class})" unless ret.kind_of? cls 1413 ret 1414 end 1415 1416 def self.coerce_to_int(obj) 1417 coerce_to(obj, Integer, :to_int) 1418 end 1419 end 1420 1421 include CoercionHelper 1422 1423 # Private CLASS 1424 1425 class Scalar < Numeric # :nodoc: 1426 include ExceptionForMatrix 1427 include CoercionHelper 1428 1429 def initialize(value) 1430 @value = value 1431 end 1432 1433 # ARITHMETIC 1434 def +(other) 1435 case other 1436 when Numeric 1437 Scalar.new(@value + other) 1438 when Vector, Matrix 1439 Scalar.Raise ErrOperationNotDefined, "+", @value.class, other.class 1440 else 1441 apply_through_coercion(other, __method__) 1442 end 1443 end 1444 1445 def -(other) 1446 case other 1447 when Numeric 1448 Scalar.new(@value - other) 1449 when Vector, Matrix 1450 Scalar.Raise ErrOperationNotDefined, "-", @value.class, other.class 1451 else 1452 apply_through_coercion(other, __method__) 1453 end 1454 end 1455 1456 def *(other) 1457 case other 1458 when Numeric 1459 Scalar.new(@value * other) 1460 when Vector, Matrix 1461 other.collect{|e| @value * e} 1462 else 1463 apply_through_coercion(other, __method__) 1464 end 1465 end 1466 1467 def / (other) 1468 case other 1469 when Numeric 1470 Scalar.new(@value / other) 1471 when Vector 1472 Scalar.Raise ErrOperationNotDefined, "/", @value.class, other.class 1473 when Matrix 1474 self * other.inverse 1475 else 1476 apply_through_coercion(other, __method__) 1477 end 1478 end 1479 1480 def ** (other) 1481 case other 1482 when Numeric 1483 Scalar.new(@value ** other) 1484 when Vector 1485 Scalar.Raise ErrOperationNotDefined, "**", @value.class, other.class 1486 when Matrix 1487 #other.powered_by(self) 1488 Scalar.Raise ErrOperationNotImplemented, "**", @value.class, other.class 1489 else 1490 apply_through_coercion(other, __method__) 1491 end 1492 end 1493 end 1494 1495end 1496 1497 1498# 1499# The +Vector+ class represents a mathematical vector, which is useful in its own right, and 1500# also constitutes a row or column of a Matrix. 1501# 1502# == Method Catalogue 1503# 1504# To create a Vector: 1505# * Vector.[](*array) 1506# * Vector.elements(array, copy = true) 1507# 1508# To access elements: 1509# * #[](i) 1510# 1511# To enumerate the elements: 1512# * #each2(v) 1513# * #collect2(v) 1514# 1515# Vector arithmetic: 1516# * #*(x) "is matrix or number" 1517# * #+(v) 1518# * #-(v) 1519# 1520# Vector functions: 1521# * #inner_product(v) 1522# * #collect 1523# * #magnitude 1524# * #map 1525# * #map2(v) 1526# * #norm 1527# * #normalize 1528# * #r 1529# * #size 1530# 1531# Conversion to other data types: 1532# * #covector 1533# * #to_a 1534# * #coerce(other) 1535# 1536# String representations: 1537# * #to_s 1538# * #inspect 1539# 1540class Vector 1541 include ExceptionForMatrix 1542 include Enumerable 1543 include Matrix::CoercionHelper 1544 extend Matrix::ConversionHelper 1545 #INSTANCE CREATION 1546 1547 private_class_method :new 1548 attr_reader :elements 1549 protected :elements 1550 1551 # 1552 # Creates a Vector from a list of elements. 1553 # Vector[7, 4, ...] 1554 # 1555 def Vector.[](*array) 1556 new convert_to_array(array, false) 1557 end 1558 1559 # 1560 # Creates a vector from an Array. The optional second argument specifies 1561 # whether the array itself or a copy is used internally. 1562 # 1563 def Vector.elements(array, copy = true) 1564 new convert_to_array(array, copy) 1565 end 1566 1567 # 1568 # Vector.new is private; use Vector[] or Vector.elements to create. 1569 # 1570 def initialize(array) 1571 # No checking is done at this point. 1572 @elements = array 1573 end 1574 1575 # ACCESSING 1576 1577 # 1578 # Returns element number +i+ (starting at zero) of the vector. 1579 # 1580 def [](i) 1581 @elements[i] 1582 end 1583 alias element [] 1584 alias component [] 1585 1586 def []=(i, v) 1587 @elements[i]= v 1588 end 1589 alias set_element []= 1590 alias set_component []= 1591 private :[]=, :set_element, :set_component 1592 1593 # 1594 # Returns the number of elements in the vector. 1595 # 1596 def size 1597 @elements.size 1598 end 1599 1600 #-- 1601 # ENUMERATIONS -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- 1602 #++ 1603 1604 # 1605 # Iterate over the elements of this vector 1606 # 1607 def each(&block) 1608 return to_enum(:each) unless block_given? 1609 @elements.each(&block) 1610 self 1611 end 1612 1613 # 1614 # Iterate over the elements of this vector and +v+ in conjunction. 1615 # 1616 def each2(v) # :yield: e1, e2 1617 raise TypeError, "Integer is not like Vector" if v.kind_of?(Integer) 1618 Vector.Raise ErrDimensionMismatch if size != v.size 1619 return to_enum(:each2, v) unless block_given? 1620 size.times do |i| 1621 yield @elements[i], v[i] 1622 end 1623 self 1624 end 1625 1626 # 1627 # Collects (as in Enumerable#collect) over the elements of this vector and +v+ 1628 # in conjunction. 1629 # 1630 def collect2(v) # :yield: e1, e2 1631 raise TypeError, "Integer is not like Vector" if v.kind_of?(Integer) 1632 Vector.Raise ErrDimensionMismatch if size != v.size 1633 return to_enum(:collect2, v) unless block_given? 1634 Array.new(size) do |i| 1635 yield @elements[i], v[i] 1636 end 1637 end 1638 1639 #-- 1640 # COMPARING -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- 1641 #++ 1642 1643 # 1644 # Returns +true+ iff the two vectors have the same elements in the same order. 1645 # 1646 def ==(other) 1647 return false unless Vector === other 1648 @elements == other.elements 1649 end 1650 1651 def eql?(other) 1652 return false unless Vector === other 1653 @elements.eql? other.elements 1654 end 1655 1656 # 1657 # Return a copy of the vector. 1658 # 1659 def clone 1660 self.class.elements(@elements) 1661 end 1662 1663 # 1664 # Return a hash-code for the vector. 1665 # 1666 def hash 1667 @elements.hash 1668 end 1669 1670 #-- 1671 # ARITHMETIC -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- 1672 #++ 1673 1674 # 1675 # Multiplies the vector by +x+, where +x+ is a number or another vector. 1676 # 1677 def *(x) 1678 case x 1679 when Numeric 1680 els = @elements.collect{|e| e * x} 1681 self.class.elements(els, false) 1682 when Matrix 1683 Matrix.column_vector(self) * x 1684 when Vector 1685 Vector.Raise ErrOperationNotDefined, "*", self.class, x.class 1686 else 1687 apply_through_coercion(x, __method__) 1688 end 1689 end 1690 1691 # 1692 # Vector addition. 1693 # 1694 def +(v) 1695 case v 1696 when Vector 1697 Vector.Raise ErrDimensionMismatch if size != v.size 1698 els = collect2(v) {|v1, v2| 1699 v1 + v2 1700 } 1701 self.class.elements(els, false) 1702 when Matrix 1703 Matrix.column_vector(self) + v 1704 else 1705 apply_through_coercion(v, __method__) 1706 end 1707 end 1708 1709 # 1710 # Vector subtraction. 1711 # 1712 def -(v) 1713 case v 1714 when Vector 1715 Vector.Raise ErrDimensionMismatch if size != v.size 1716 els = collect2(v) {|v1, v2| 1717 v1 - v2 1718 } 1719 self.class.elements(els, false) 1720 when Matrix 1721 Matrix.column_vector(self) - v 1722 else 1723 apply_through_coercion(v, __method__) 1724 end 1725 end 1726 1727 # 1728 # Vector division. 1729 # 1730 def /(x) 1731 case x 1732 when Numeric 1733 els = @elements.collect{|e| e / x} 1734 self.class.elements(els, false) 1735 when Matrix, Vector 1736 Vector.Raise ErrOperationNotDefined, "/", self.class, x.class 1737 else 1738 apply_through_coercion(x, __method__) 1739 end 1740 end 1741 1742 #-- 1743 # VECTOR FUNCTIONS -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- 1744 #++ 1745 1746 # 1747 # Returns the inner product of this vector with the other. 1748 # Vector[4,7].inner_product Vector[10,1] => 47 1749 # 1750 def inner_product(v) 1751 Vector.Raise ErrDimensionMismatch if size != v.size 1752 1753 p = 0 1754 each2(v) {|v1, v2| 1755 p += v1 * v2.conj 1756 } 1757 p 1758 end 1759 1760 # 1761 # Like Array#collect. 1762 # 1763 def collect(&block) # :yield: e 1764 return to_enum(:collect) unless block_given? 1765 els = @elements.collect(&block) 1766 self.class.elements(els, false) 1767 end 1768 alias map collect 1769 1770 # 1771 # Returns the modulus (Pythagorean distance) of the vector. 1772 # Vector[5,8,2].r => 9.643650761 1773 # 1774 def magnitude 1775 Math.sqrt(@elements.inject(0) {|v, e| v + e.abs2}) 1776 end 1777 alias r magnitude 1778 alias norm magnitude 1779 1780 # 1781 # Like Vector#collect2, but returns a Vector instead of an Array. 1782 # 1783 def map2(v, &block) # :yield: e1, e2 1784 return to_enum(:map2, v) unless block_given? 1785 els = collect2(v, &block) 1786 self.class.elements(els, false) 1787 end 1788 1789 class ZeroVectorError < StandardError 1790 end 1791 # 1792 # Returns a new vector with the same direction but with norm 1. 1793 # v = Vector[5,8,2].normalize 1794 # # => Vector[0.5184758473652127, 0.8295613557843402, 0.20739033894608505] 1795 # v.norm => 1.0 1796 # 1797 def normalize 1798 n = magnitude 1799 raise ZeroVectorError, "Zero vectors can not be normalized" if n == 0 1800 self / n 1801 end 1802 1803 #-- 1804 # CONVERTING 1805 #++ 1806 1807 # 1808 # Creates a single-row matrix from this vector. 1809 # 1810 def covector 1811 Matrix.row_vector(self) 1812 end 1813 1814 # 1815 # Returns the elements of the vector in an array. 1816 # 1817 def to_a 1818 @elements.dup 1819 end 1820 1821 def elements_to_f 1822 warn "#{caller(1)[0]}: warning: Vector#elements_to_f is deprecated" 1823 map(&:to_f) 1824 end 1825 1826 def elements_to_i 1827 warn "#{caller(1)[0]}: warning: Vector#elements_to_i is deprecated" 1828 map(&:to_i) 1829 end 1830 1831 def elements_to_r 1832 warn "#{caller(1)[0]}: warning: Vector#elements_to_r is deprecated" 1833 map(&:to_r) 1834 end 1835 1836 # 1837 # The coerce method provides support for Ruby type coercion. 1838 # This coercion mechanism is used by Ruby to handle mixed-type 1839 # numeric operations: it is intended to find a compatible common 1840 # type between the two operands of the operator. 1841 # See also Numeric#coerce. 1842 # 1843 def coerce(other) 1844 case other 1845 when Numeric 1846 return Matrix::Scalar.new(other), self 1847 else 1848 raise TypeError, "#{self.class} can't be coerced into #{other.class}" 1849 end 1850 end 1851 1852 #-- 1853 # PRINTING -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- 1854 #++ 1855 1856 # 1857 # Overrides Object#to_s 1858 # 1859 def to_s 1860 "Vector[" + @elements.join(", ") + "]" 1861 end 1862 1863 # 1864 # Overrides Object#inspect 1865 # 1866 def inspect 1867 "Vector" + @elements.inspect 1868 end 1869end 1870