← Back to branch summary

~ubuntu-branches/ubuntu/quantal/ruby1.9.1/quantal

« back to all changes in this revision

Viewing changes to lib/matrix.rb

  • Committer: Bazaar Package Importer
  • Author(s): Lucas Nussbaum
  • Date: 2011-09-24 19:16:17 UTC
  • mfrom: (1.1.8 upstream) (13.1.7 experimental)
  • Revision ID: james.westby@ubuntu.com-20110924191617-o1qz4rcmqjot8zuy
Tags: 1.9.3~rc1-1
http://bugs.debian.org/639959
* New upstream release: 1.9.3 RC1.
  + Includes load.c fixes. Closes: #639959.
* Upload to unstable.

Show diffs side-by-side

added added

removed removed

Lines of Context:
lib/matrix.rb
56
56
# * <tt> #map                           </tt>
57
57
# * <tt> #each                          </tt>
58
58
# * <tt> #each_with_index               </tt>
 
59
# * <tt> #find_index                    </tt>
59
60
# * <tt> #minor(*param)                 </tt>
60
61
#
61
62
# Properties of a matrix:
 
63
# * <tt> #diagonal?                     </tt>
62
64
# * <tt> #empty?                        </tt>
 
65
# * <tt> #hermitian?                    </tt>
 
66
# * <tt> #lower_triangular?             </tt>
 
67
# * <tt> #normal?                       </tt>
 
68
# * <tt> #orthogonal?                   </tt>
 
69
# * <tt> #permutation?                  </tt>
63
70
# * <tt> #real?                         </tt>
64
71
# * <tt> #regular?                      </tt>
65
72
# * <tt> #singular?                     </tt>
66
73
# * <tt> #square?                       </tt>
 
74
# * <tt> #symmetric?                    </tt>
 
75
# * <tt> #unitary?                      </tt>
 
76
# * <tt> #upper_triangular?             </tt>
 
77
# * <tt> #zero?                         </tt>
67
78
#
68
79
# Matrix arithmetic:
69
80
# * <tt>  *(m)                          </tt>
78
89
# * <tt> #determinant                   </tt>
79
90
# * <tt> #det                           </tt>
80
91
# * <tt> #rank                          </tt>
 
92
# * <tt> #round                         </tt>
81
93
# * <tt> #trace                         </tt>
82
94
# * <tt> #tr                            </tt>
83
95
# * <tt> #transpose                     </tt>
84
96
# * <tt> #t                             </tt>
85
97
#
 
98
# Matrix decompositions:
 
99
# * <tt> #eigen                         </tt>
 
100
# * <tt> #eigensystem                   </tt>
 
101
# * <tt> #lup                           </tt>
 
102
# * <tt> #lup_decomposition             </tt>
 
103
#
86
104
# Complex arithmetic:
87
105
# * <tt> conj                           </tt>
88
106
# * <tt> conjugate                      </tt>
105
123
class Matrix
106
124
  include Enumerable
107
125
  include ExceptionForMatrix
 
126
  autoload :EigenvalueDecomposition, "matrix/eigenvalue_decomposition"
 
127
  autoload :LUPDecomposition, "matrix/lup_decomposition"
108
128
 
109
129
  # instance creations
110
130
  private_class_method :new
218
238
  end
219
239
 
220
240
  #
221
 
  # Creates an +n+ by +n+ zero matrix.
 
241
  # Creates a zero matrix.
222
242
  #   Matrix.zero(2)
223
243
  #     => 0 0
224
244
  #        0 0
225
245
  #
226
 
  def Matrix.zero(n)
227
 
    Matrix.scalar(n, 0)
 
246
  def Matrix.zero(row_size, column_size = row_size)
 
247
    rows = Array.new(row_size){Array.new(column_size, 0)}
 
248
    new rows, column_size
228
249
  end
229
250
 
230
251
  #
365
386
 
366
387
  #
367
388
  # Yields all elements of the matrix, starting with those of the first row,
368
 
  # or returns an Enumerator is no block given
 
389
  # or returns an Enumerator is no block given.
 
390
  # Elements can be restricted by passing an argument:
 
391
  # * :all (default): yields all elements
 
392
  # * :diagonal: yields only elements on the diagonal
 
393
  # * :off_diagonal: yields all elements except on the diagonal
 
394
  # * :lower: yields only elements on or below the diagonal
 
395
  # * :strict_lower: yields only elements below the diagonal
 
396
  # * :strict_upper: yields only elements above the diagonal
 
397
  # * :upper: yields only elements on or above the diagonal
 
398
  #
369
399
  #   Matrix[ [1,2], [3,4] ].each { |e| puts e }
370
400
  #     # => prints the numbers 1 to 4
 
401
  #   Matrix[ [1,2], [3,4] ].each(:strict_lower).to_a # => [3]
371
402
  #
372
 
  def each(&block) # :yield: e
373
 
    return to_enum(:each) unless block_given?
374
 
    @rows.each do |row|
375
 
      row.each(&block)
 
403
  def each(which = :all) # :yield: e
 
404
    return to_enum :each, which unless block_given?
 
405
    last = column_size - 1
 
406
    case which
 
407
    when :all
 
408
      block = Proc.new
 
409
      @rows.each do |row|
 
410
        row.each(&block)
 
411
      end
 
412
    when :diagonal
 
413
      @rows.each_with_index do |row, row_index|
 
414
        yield row.fetch(row_index){return self}
 
415
      end
 
416
    when :off_diagonal
 
417
      @rows.each_with_index do |row, row_index|
 
418
        column_size.times do |col_index|
 
419
          yield row[col_index] unless row_index == col_index
 
420
        end
 
421
      end
 
422
    when :lower
 
423
      @rows.each_with_index do |row, row_index|
 
424
        0.upto([row_index, last].min) do |col_index|
 
425
          yield row[col_index]
 
426
        end
 
427
      end
 
428
    when :strict_lower
 
429
      @rows.each_with_index do |row, row_index|
 
430
        [row_index, column_size].min.times do |col_index|
 
431
          yield row[col_index]
 
432
        end
 
433
      end
 
434
    when :strict_upper
 
435
      @rows.each_with_index do |row, row_index|
 
436
        (row_index+1).upto(last) do |col_index|
 
437
          yield row[col_index]
 
438
        end
 
439
      end
 
440
    when :upper
 
441
      @rows.each_with_index do |row, row_index|
 
442
        row_index.upto(last) do |col_index|
 
443
          yield row[col_index]
 
444
        end
 
445
      end
 
446
    else
 
447
      Matrix.Raise ArgumentError, "expected #{which.inspect} to be one of :all, :diagonal, :off_diagonal, :lower, :strict_lower, :strict_upper or :upper"
376
448
    end
377
449
    self
378
450
  end
379
451
 
380
452
  #
381
 
  # Yields all elements of the matrix, starting with those of the first row,
382
 
  # along with the row index and column index,
383
 
  # or returns an Enumerator is no block given
 
453
  # Same as #each, but the row index and column index in addition to the element
 
454
  #
384
455
  #   Matrix[ [1,2], [3,4] ].each_with_index do |e, row, col|
385
456
  #     puts "#{e} at #{row}, #{col}"
386
457
  #   end
387
 
  #     # => 1 at 0, 0
388
 
  #     # => 2 at 0, 1
389
 
  #     # => 3 at 1, 0
390
 
  #     # => 4 at 1, 1
 
458
  #     # => Prints:
 
459
  #     #    1 at 0, 0
 
460
  #     #    2 at 0, 1
 
461
  #     #    3 at 1, 0
 
462
  #     #    4 at 1, 1
391
463
  #
392
 
  def each_with_index(&block) # :yield: e, row, column
393
 
    return to_enum(:each_with_index) unless block_given?
394
 
    @rows.each_with_index do |row, row_index|
395
 
      row.each_with_index do |e, col_index|
396
 
        yield e, row_index, col_index
397
 
      end
 
464
  def each_with_index(which = :all) # :yield: e, row, column
 
465
    return to_enum :each_with_index, which unless block_given?
 
466
    last = column_size - 1
 
467
    case which
 
468
    when :all
 
469
      @rows.each_with_index do |row, row_index|
 
470
        row.each_with_index do |e, col_index|
 
471
          yield e, row_index, col_index
 
472
        end
 
473
      end
 
474
    when :diagonal
 
475
      @rows.each_with_index do |row, row_index|
 
476
        yield row.fetch(row_index){return self}, row_index, row_index
 
477
      end
 
478
    when :off_diagonal
 
479
      @rows.each_with_index do |row, row_index|
 
480
        column_size.times do |col_index|
 
481
          yield row[col_index], row_index, col_index unless row_index == col_index
 
482
        end
 
483
      end
 
484
    when :lower
 
485
      @rows.each_with_index do |row, row_index|
 
486
        0.upto([row_index, last].min) do |col_index|
 
487
          yield row[col_index], row_index, col_index
 
488
        end
 
489
      end
 
490
    when :strict_lower
 
491
      @rows.each_with_index do |row, row_index|
 
492
        [row_index, column_size].min.times do |col_index|
 
493
          yield row[col_index], row_index, col_index
 
494
        end
 
495
      end
 
496
    when :strict_upper
 
497
      @rows.each_with_index do |row, row_index|
 
498
        (row_index+1).upto(last) do |col_index|
 
499
          yield row[col_index], row_index, col_index
 
500
        end
 
501
      end
 
502
    when :upper
 
503
      @rows.each_with_index do |row, row_index|
 
504
        row_index.upto(last) do |col_index|
 
505
          yield row[col_index], row_index, col_index
 
506
        end
 
507
      end
 
508
    else
 
509
      Matrix.Raise ArgumentError, "expected #{which.inspect} to be one of :all, :diagonal, :off_diagonal, :lower, :strict_lower, :strict_upper or :upper"
398
510
    end
399
511
    self
400
512
  end
401
513
 
 
514
  SELECTORS = {all: true, diagonal: true, off_diagonal: true, lower: true, strict_lower: true, strict_upper: true, upper: true}.freeze
 
515
  #
 
516
  # :call-seq:
 
517
  #   index(value, selector = :all) -> [row, column]
 
518
  #   index(selector = :all){ block } -> [row, column]
 
519
  #   index(selector = :all) -> an_enumerator
 
520
  #
 
521
  # The index method is specialized to return the index as [row, column]
 
522
  # It also accepts an optional +selector+ argument, see #each for details.
 
523
  #
 
524
  #   Matrix[ [1,2], [3,4] ].index(&:even?) # => [0, 1]
 
525
  #   Matrix[ [1,1], [1,1] ].index(1, :strict_lower) # => [1, 0]
 
526
  #
 
527
  def index(*args)
 
528
    raise ArgumentError, "wrong number of arguments(#{args.size} for 0-2)" if args.size > 2
 
529
    which = (args.size == 2 || SELECTORS.include?(args.last)) ? args.pop : :all
 
530
    return to_enum :find_index, which, *args unless block_given? || args.size == 1
 
531
    if args.size == 1
 
532
      value = args.first
 
533
      each_with_index(which) do |e, row_index, col_index|
 
534
        return row_index, col_index if e == value
 
535
      end
 
536
    else
 
537
      each_with_index(which) do |e, row_index, col_index|
 
538
        return row_index, col_index if yield e
 
539
      end
 
540
    end
 
541
    nil
 
542
  end
 
543
  alias_method :find_index, :index
402
544
  #
403
545
  # Returns a section of the matrix.  The parameters are either:
404
546
  # *  start_row, nrows, start_col, ncols; OR
450
592
  #++
451
593
 
452
594
  #
 
595
  # Returns +true+ is this is a diagonal matrix.
 
596
  # Raises an error if matrix is not square.
 
597
  #
 
598
  def diagonal?
 
599
    Matrix.Raise ErrDimensionMismatch unless square?
 
600
    each(:off_diagonal).all?(&:zero?)
 
601
  end
 
602
 
 
603
  #
453
604
  # Returns +true+ if this is an empty matrix, i.e. if the number of rows
454
605
  # or the number of columns is 0.
455
606
  #
458
609
  end
459
610
 
460
611
  #
 
612
  # Returns +true+ is this is an hermitian matrix.
 
613
  # Raises an error if matrix is not square.
 
614
  #
 
615
  def hermitian?
 
616
    Matrix.Raise ErrDimensionMismatch unless square?
 
617
    each_with_index(:strict_upper).all? do |e, row, col|
 
618
      e == rows[col][row].conj
 
619
    end
 
620
  end
 
621
 
 
622
  #
 
623
  # Returns +true+ is this is a lower triangular matrix.
 
624
  #
 
625
  def lower_triangular?
 
626
    each(:strict_upper).all?(&:zero?)
 
627
  end
 
628
 
 
629
  #
 
630
  # Returns +true+ is this is a normal matrix.
 
631
  # Raises an error if matrix is not square.
 
632
  #
 
633
  def normal?
 
634
    Matrix.Raise ErrDimensionMismatch unless square?
 
635
    rows.each_with_index do |row_i, i|
 
636
      rows.each_with_index do |row_j, j|
 
637
        s = 0
 
638
        rows.each_with_index do |row_k, k|
 
639
          s += row_i[k] * row_j[k].conj - row_k[i].conj * row_k[j]
 
640
        end
 
641
        return false unless s == 0
 
642
      end
 
643
    end
 
644
    true
 
645
  end
 
646
 
 
647
  #
 
648
  # Returns +true+ is this is an orthogonal matrix
 
649
  # Raises an error if matrix is not square.
 
650
  #
 
651
  def orthogonal?
 
652
    Matrix.Raise ErrDimensionMismatch unless square?
 
653
    rows.each_with_index do |row, i|
 
654
      column_size.times do |j|
 
655
        s = 0
 
656
        row_size.times do |k|
 
657
          s += row[k] * rows[k][j]
 
658
        end
 
659
        return false unless s == (i == j ? 1 : 0)
 
660
      end
 
661
    end
 
662
    true
 
663
  end
 
664
 
 
665
  #
 
666
  # Returns +true+ is this is a permutation matrix
 
667
  # Raises an error if matrix is not square.
 
668
  #
 
669
  def permutation?
 
670
    Matrix.Raise ErrDimensionMismatch unless square?
 
671
    cols = Array.new(column_size)
 
672
    rows.each_with_index do |row, i|
 
673
      found = false
 
674
      row.each_with_index do |e, j|
 
675
        if e == 1
 
676
          return false if found || cols[j]
 
677
          found = cols[j] = true
 
678
        elsif e != 0
 
679
          return false
 
680
        end
 
681
      end
 
682
      return false unless found
 
683
    end
 
684
    true
 
685
  end
 
686
 
 
687
  #
461
688
  # Returns +true+ if all entries of the matrix are real.
462
689
  #
463
690
  def real?
485
712
    column_size == row_size
486
713
  end
487
714
 
 
715
  #
 
716
  # Returns +true+ is this is a symmetric matrix.
 
717
  # Raises an error if matrix is not square.
 
718
  #
 
719
  def symmetric?
 
720
    Matrix.Raise ErrDimensionMismatch unless square?
 
721
    each_with_index(:strict_upper).all? do |e, row, col|
 
722
      e == rows[col][row]
 
723
    end
 
724
  end
 
725
 
 
726
  #
 
727
  # Returns +true+ is this is a unitary matrix
 
728
  # Raises an error if matrix is not square.
 
729
  #
 
730
  def unitary?
 
731
    Matrix.Raise ErrDimensionMismatch unless square?
 
732
    rows.each_with_index do |row, i|
 
733
      column_size.times do |j|
 
734
        s = 0
 
735
        row_size.times do |k|
 
736
          s += row[k].conj * rows[k][j]
 
737
        end
 
738
        return false unless s == (i == j ? 1 : 0)
 
739
      end
 
740
    end
 
741
    true
 
742
  end
 
743
 
 
744
  #
 
745
  # Returns +true+ is this is an upper triangular matrix.
 
746
  #
 
747
  def upper_triangular?
 
748
    each(:strict_lower).all?(&:zero?)
 
749
  end
 
750
 
 
751
  #
 
752
  # Returns +true+ is this is a matrix with only zero elements
 
753
  #
 
754
  def zero?
 
755
    all?(&:zero?)
 
756
  end
 
757
 
488
758
  #--
489
759
  # OBJECT METHODS -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
490
760
  #++
689
959
  private :inverse_from
690
960
 
691
961
  #
692
 
  # Matrix exponentiation.  Currently implemented for integer powers only.
 
962
  # Matrix exponentiation.
693
963
  # Equivalent to multiplying the matrix by itself N times.
 
964
  # Non integer exponents will be handled by diagonalizing the matrix.
 
965
  #
694
966
  #   Matrix[[7,6], [3,9]] ** 2
695
967
  #     => 67 96
696
968
  #        48 99
710
982
        return z if (other >>= 1).zero?
711
983
        x *= x
712
984
      end
713
 
    when Float, Rational
714
 
      Matrix.Raise ErrOperationNotImplemented, "**", self.class, other.class
 
985
    when Numeric
 
986
      v, d, v_inv = eigensystem
 
987
      v * Matrix.diagonal(*d.each(:diagonal).map{|e| e ** other}) * v_inv
715
988
    else
716
989
      Matrix.Raise ErrOperationNotDefined, "**", self.class, other.class
717
990
    end
834
1107
    a = to_a
835
1108
    last_column = column_size - 1
836
1109
    last_row = row_size - 1
837
 
    rank = 0
838
1110
    pivot_row = 0
839
1111
    previous_pivot = 1
840
1112
    0.upto(last_column) do |k|
865
1137
    rank
866
1138
  end
867
1139
 
 
1140
  # Returns a matrix with entries rounded to the given precision
 
1141
  # (see Float#round)
 
1142
  #
 
1143
  def round(ndigits=0)
 
1144
    map{|e| e.round(ndigits)}
 
1145
  end
868
1146
 
869
1147
  #
870
1148
  # Returns the trace (sum of diagonal elements) of the matrix.
896
1174
  alias t transpose
897
1175
 
898
1176
  #--
 
1177
  # DECOMPOSITIONS -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
 
1178
  #++
 
1179
 
 
1180
  #
 
1181
  # Returns the Eigensystem of the matrix; see +EigenvalueDecomposition+.
 
1182
  #   m = Matrix[[1, 2], [3, 4]]
 
1183
  #   v, d, v_inv = m.eigensystem
 
1184
  #   d.diagonal? # => true
 
1185
  #   v.inv == v_inv # => true
 
1186
  #   (v * d * v_inv).round(5) == m # => true
 
1187
  #
 
1188
  def eigensystem
 
1189
    EigenvalueDecomposition.new(self)
 
1190
  end
 
1191
  alias eigen eigensystem
 
1192
 
 
1193
  #
 
1194
  # Returns the LUP decomposition of the matrix; see +LUPDecomposition+.
 
1195
  #   a = Matrix[[1, 2], [3, 4]]
 
1196
  #   l, u, p = a.lup
 
1197
  #   l.lower_triangular? # => true
 
1198
  #   u.upper_triangular? # => true
 
1199
  #   p.permutation?      # => true
 
1200
  #   l * u == a * p      # => true
 
1201
  #   a.lup.solve([2, 5]) # => Vector[(1/1), (1/2)]
 
1202
  #
 
1203
  def lup
 
1204
    LUPDecomposition.new(self)
 
1205
  end
 
1206
  alias lup_decomposition lup
 
1207
 
 
1208
  #--
899
1209
  # COMPLEX ARITHMETIC -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
900
1210
  #++
901
1211
 
1207
1517
# Vector functions:
1208
1518
# * <tt> #inner_product(v)                    </tt>
1209
1519
# * <tt> #collect                             </tt>
 
1520
# * <tt> #magnitude                           </tt>
1210
1521
# * <tt> #map                                 </tt>
1211
1522
# * <tt> #map2(v)                             </tt>
 
1523
# * <tt> #norm                                </tt>
 
1524
# * <tt> #normalize                           </tt>
1212
1525
# * <tt> #r                                   </tt>
1213
1526
# * <tt> #size                                </tt>
1214
1527
#
1237
1550
  #   Vector[7, 4, ...]
1238
1551
  #
1239
1552
  def Vector.[](*array)
1240
 
    new convert_to_array(array, copy = false)
 
1553
    new convert_to_array(array, false)
1241
1554
  end
1242
1555
 
1243
1556
  #
1452
1765
  alias map collect
1453
1766
 
1454
1767
  #
 
1768
  # Returns the modulus (Pythagorean distance) of the vector.
 
1769
  #   Vector[5,8,2].r => 9.643650761
 
1770
  #
 
1771
  def magnitude
 
1772
    Math.sqrt(@elements.inject(0) {|v, e| v + e*e})
 
1773
  end
 
1774
  alias r magnitude
 
1775
  alias norm magnitude
 
1776
 
 
1777
  #
1455
1778
  # Like Vector#collect2, but returns a Vector instead of an Array.
1456
1779
  #
1457
1780
  def map2(v, &block) # :yield: e1, e2
1460
1783
    Vector.elements(els, false)
1461
1784
  end
1462
1785
 
1463
 
  #
1464
 
  # Returns the modulus (Pythagorean distance) of the vector.
1465
 
  #   Vector[5,8,2].r => 9.643650761
1466
 
  #
1467
 
  def r
1468
 
    Math.sqrt(@elements.inject(0) {|v, e| v + e*e})
 
1786
  class ZeroVectorError < StandardError
 
1787
  end
 
1788
  #
 
1789
  # Returns a new vector with the same direction but with norm 1.
 
1790
  #   v = Vector[5,8,2].normalize
 
1791
  #   # => Vector[0.5184758473652127, 0.8295613557843402, 0.20739033894608505]
 
1792
  #   v.norm => 1.0
 
1793
  #
 
1794
  def normalize
 
1795
    n = magnitude
 
1796
    raise ZeroVectorError, "Zero vectors can not be normalized" if n == 0
 
1797
    self / n
1469
1798
  end
1470
1799
 
1471
1800
  #--
1532
1861
  # Overrides Object#inspect
1533
1862
  #
1534
1863
  def inspect
1535
 
    str = "Vector"+@elements.inspect
 
1864
    "Vector" + @elements.inspect
1536
1865
  end
1537
1866
end