lib/matrix.rb
DEFINITIONS
This source file includes following functions.
1 #!/usr/local/bin/ruby
2 #
3 # matrix.rb -
4 # $Release Version: 1.0$
5 # $Revision: 1.11 $
6 # $Date: 1999/10/06 11:01:53 $
7 # Original Version from Smalltalk-80 version
8 # on July 23, 1985 at 8:37:17 am
9 # by Keiju ISHITSUKA
10 #
11 # --
12 #
13 # Matrix[[1,2,3],
14 # :
15 # [3,4,5]]
16 # Matrix[row0,
17 # row1,
18 # :
19 # rown]
20 #
21 #
22 # module ExceptionForMatrix::
23 # Exceptions:
24 # ErrDimensionMismatch
25 # number of column/row do not match
26 # ErrNotRegular
27 # not a regular matrix
28 # ErrOperationNotDefined
29 # specified operator is not defined (yet)
30 #
31 # class Matrix
32 # include ExceptionForMatrix
33 #
34 # Methods:
35 # class methods:
36 # Matrix.[](*rows)
37 # creates a matrix where `rows' indicates rows.
38 # `rows' is an array of arrays,
39 # e.g, Matrix[[11, 12], [21, 22]]
40 # Matrix.rows(rows, copy = true)
41 # creates a matrix where `rows' indicates rows.
42 # if optional argument `copy' is false, use the array as
43 # internal structure of the metrix without copying.
44 # Matrix.columns(columns)
45 # creates a new matrix using `columns` as set of colums vectors.
46 # Matrix.diagonal(*values)
47 # creates a matrix where `columns' indicates columns.
48 # Matrix.scalar(n, value)
49 # creates a diagonal matrix such that the diagal compornents is
50 # given by `values'.
51 # Matrix.scalar(n, value)
52 # creates an n-by-n scalar matrix such that the diagal compornent is
53 # given by `value'.
54 # Matrix.identity(n)
55 # Matrix.unit(n)
56 # Matrix.I(n)
57 # creates an n-by-n unit matrix.
58 # Matrix.zero(n)
59 # creates an n-by-n zero matrix.
60 # Matrix.row_vector(row)
61 # creates a 1-by-n matrix such the row vector is `row'.
62 # `row' is specifed as a Vector or an Array.
63 # Matrix.column_vector(column)
64 # creates a 1-by-n matrix such that column vector is `column'.
65 # `column' is specifed as a Vector or an Array.
66 # accessing:
67 # [](i, j)
68 # returns (i,j) compornent
69 # row_size
70 # returns the number of rows
71 # column_size
72 # returns the number of columns
73 # row(i)
74 # returns the i-th row vector.
75 # when the block is supplied for the method, the block is iterated
76 # over all row vectors.
77 # column(j)
78 # returns the jth column vector.
79 # when the block is supplied for the method, the block is iterated
80 # over all column vectors.
81 # collect
82 # map
83 # creates a matrix which is the result of iteration of given
84 # block over all compornents.
85 # minor(*param)
86 # returns sub matrix. parameter is specified as the following:
87 # 1. from_row, row_size, from_col, size_col
88 # 2. from_row..to_row, from_col..to_col
89 # TESTING:
90 # regular?
91 # Is regular?
92 # singular?
93 # Is singular? i.e. Is non-regular?
94 # square?
95 # Is square?
96 # ARITHMETIC:
97 # *(m)
98 # times
99 # +(m)
100 # plus
101 # -(m)
102 # minus
103 # /(m)
104 # self * m.inv
105 # inverse
106 # inv
107 # inverse
108 # **
109 # power
110 # Matrix functions:
111 # determinant
112 # det
113 # returns the determinant
114 # rank
115 # returns the rank
116 # trace
117 # tr
118 # returns the trace
119 # transpose
120 # t
121 # returns the transposed
122 # CONVERTING:
123 # coerce(other)
124 # row_vectors
125 # array of row vectors
126 # column_vectors
127 # array of column vectors
128 # to_a
129 # converts each element to Array
130 # to_f
131 # converts each element to Float
132 # to_i
133 # converts each element to Integer
134 # to_r
135 # converts each element to Rational
136 # PRINTING:
137 # to_s
138 # returns string representation
139 # inspect
140 #
141 # class Vector
142 # include ExceptionForMatrix
143 #
144 # INSTANCE CREATION:
145 # Vector.[](*array)
146 # Vector.elements(array, copy = true)
147 # ACCSESSING:
148 # [](i)
149 # size
150 # ENUMRATIONS:
151 # each2(v)
152 # collect2(v)
153 # ARITHMETIC:
154 # *(x) "is matrix or number"
155 # +(v)
156 # -(v)
157 # VECTOR FUNCTIONS:
158 # inner_product(v)
159 # collect
160 # map
161 # map2(v)
162 # r
163 # CONVERTING:
164 # covector
165 # to_a
166 # to_f
167 # to_i
168 # to_r
169 # coerce(other)
170 # PRINTING:
171 # to_s
172 # inspect
173
174 require "e2mmap.rb"
175
176 module ExceptionForMatrix
177 extend Exception2MessageMapper
178 def_e2message(TypeError, "wrong argument type %s (expected %s)")
179 def_e2message(ArgumentError, "Wrong # of arguments(%d for %d)")
180
181 def_exception("ErrDimensionMismatch", "\#{self.name} dimension mismatch")
182 def_exception("ErrNotRegular", "Not Regular Matrix")
183 def_exception("ErrOperationNotDefined", "This operation(%s) can\\'t defined")
184 end
185
186 class Matrix
187 @RCS_ID='-$Id: matrix.rb,v 1.11 1999/10/06 11:01:53 keiju Exp keiju $-'
188
189 # extend Exception2MessageMapper
190 include ExceptionForMatrix
191
192 # instance creations
193 private_class_method :new
194
195 def Matrix.[](*rows)
196 new(:init_rows, rows, false)
197 end
198
199 def Matrix.rows(rows, copy = true)
200 new(:init_rows, rows, copy)
201 end
202
203 def Matrix.columns(columns)
204 rows = (0 .. columns[0].size - 1).collect {
205 |i|
206 (0 .. columns.size - 1).collect {
207 |j|
208 columns[j][i]
209 }
210 }
211 Matrix.rows(rows, false)
212 end
213
214 def Matrix.diagonal(*values)
215 size = values.size
216 rows = (0 .. size - 1).collect {
217 |j|
218 row = Array.new(size).fill(0, 0, size)
219 row[j] = values[j]
220 row
221 }
222 rows(rows, false)
223 end
224
225 def Matrix.scalar(n, value)
226 Matrix.diagonal(*Array.new(n).fill(value, 0, n))
227 end
228
229 def Matrix.identity(n)
230 Matrix.scalar(n, 1)
231 end
232 class << Matrix
233 alias unit identity
234 alias I identity
235 end
236
237 def Matrix.zero(n)
238 Matrix.scalar(n, 0)
239 end
240
241 def Matrix.row_vector(row)
242 case row
243 when Vector
244 Matrix.rows([row.to_a], false)
245 when Array
246 Matrix.rows([row.dup], false)
247 else
248 Matrix.row([[row]], false)
249 end
250 end
251
252 def Matrix.column_vector(column)
253 case column
254 when Vector
255 Matrix.columns([column.to_a])
256 when Array
257 Matrix.columns([column])
258 else
259 Matrix.columns([[column]])
260 end
261 end
262
263 # initializing
264 def initialize(init_method, *argv)
265 self.send(init_method, *argv)
266 end
267
268 def init_rows(rows, copy)
269 if copy
270 @rows = rows.collect{|row| row.dup}
271 else
272 @rows = rows
273 end
274 self
275 end
276 private :init_rows
277
278 #accessing
279 def [](i, j)
280 @rows[i][j]
281 end
282
283 def row_size
284 @rows.size
285 end
286
287 def column_size
288 @rows[0].size
289 end
290
291 def row(i)
292 if block_given?
293 for e in @rows[i]
294 yield e
295
296 end
297 else
298 Vector.elements(@rows[i])
299 end
300 end
301
302 def column(j)
303 if block_given?
304 0.upto(row_size - 1) do
305 |i|
306 yield @rows[i][j]
307 end
308 else
309 col = (0 .. row_size - 1).collect {
310 |i|
311 @rows[i][j]
312 }
313 Vector.elements(col, false)
314 end
315 end
316
317 def collect
318 rows = @rows.collect{|row| row.collect{|e| yield e}}
319 Matrix.rows(rows, false)
320 end
321 alias map collect
322
323 #
324 # param: (from_row, row_size, from_col, size_col)
325 # (from_row..to_row, from_col..to_col)
326 #
327 def minor(*param)
328 case param.size
329 when 2
330 from_row = param[0].first
331 size_row = param[0].size
332 from_col = param[1].first
333 size_col = param[1].size
334 when 4
335 from_row = param[0]
336 size_row = param[1]
337 from_col = param[2]
338 size_col = param[3]
339 else
340 Matrix.Raise ArgumentError, param.inspect
341 end
342
343 rows = @rows[from_row, size_row].collect{
344 |row|
345 row[from_col, size_col]
346 }
347 Matrix.rows(rows, false)
348 end
349
350 # TESTING
351 def regular?
352 square? and rank == column_size
353 end
354
355 def singular?
356 not regular?
357 end
358
359 def square?
360 column_size == row_size
361 end
362
363 # COMPARING
364 def ==(other)
365 return false unless Matrix === other
366
367 other.compare_by_row_vectors(@rows)
368 end
369 alias eql? ==
370
371 def compare_by_row_vectors(rows)
372 return false unless @rows.size == rows.size
373
374 0.upto(@rows.size - 1) do
375 |i|
376 return false unless @rows[i] == rows[i]
377 end
378 true
379 end
380
381 def clone
382 Matrix.rows(@rows)
383 end
384
385 def hash
386 value = 0
387 for row in @rows
388 for e in row
389 value ^= e.hash
390 end
391 end
392 return value
393 end
394
395 # ARITHMETIC
396
397 def *(m) # m is matrix or vector or number
398 case(m)
399 when Numeric
400 rows = @rows.collect {
401 |row|
402 row.collect {
403 |e|
404 e * m
405 }
406 }
407 return Matrix.rows(rows, false)
408 when Vector
409 m = Matrix.column_vector(m)
410 r = self * m
411 return r.column(0)
412 when Matrix
413 Matrix.Raise ErrDimensionMismatch if column_size != m.row_size
414
415 rows = (0 .. row_size - 1).collect {
416 |i|
417 (0 .. m.column_size - 1).collect {
418 |j|
419 vij = 0
420 0.upto(column_size - 1) do
421 |k|
422 vij += self[i, k] * m[k, j]
423 end
424 vij
425 }
426 }
427 return Matrix.rows(rows, false)
428 else
429 x, y = m.coerce(self)
430 return x * y
431 end
432 end
433
434 def +(m)
435 case m
436 when Numeric
437 Matrix.Raise ErrOperationNotDefined, "+"
438 when Vector
439 m = Matrix.column_vector(m)
440 when Matrix
441 else
442 x, y = m.coerce(self)
443 return x + y
444 end
445
446 Matrix.Raise ErrDimensionMismatch unless row_size == m.row_size and column_size == m.column_size
447
448 rows = (0 .. row_size - 1).collect {
449 |i|
450 (0 .. column_size - 1).collect {
451 |j|
452 self[i, j] + m[i, j]
453 }
454 }
455 Matrix.rows(rows, false)
456 end
457
458 def -(m)
459 case m
460 when Numeric
461 Matrix.Raise ErrOperationNotDefined, "-"
462 when Vector
463 m = Matrix.column_vector(m)
464 when Matrix
465 else
466 x, y = m.coerce(self)
467 return x - y
468 end
469
470 Matrix.Raise ErrDimensionMismatch unless row_size == m.row_size and column_size == m.column_size
471
472 rows = (0 .. row_size - 1).collect {
473 |i|
474 (0 .. column_size - 1).collect {
475 |j|
476 self[i, j] - m[i, j]
477 }
478 }
479 Matrix.rows(rows, false)
480 end
481
482 def /(other)
483 case other
484 when Numeric
485 rows = @rows.collect {
486 |row|
487 row.collect {
488 |e|
489 e / other
490 }
491 }
492 return Matrix.rows(rows, false)
493 when Matrix
494 return self * other.inverse
495 else
496 x, y = other.coerce(self)
497 rerurn x / y
498 end
499 end
500
501 def inverse
502 Matrix.Raise ErrDimensionMismatch unless square?
503 Matrix.I(row_size).inverse_from(self)
504 end
505 alias inv inverse
506
507 def inverse_from(src)
508 size = row_size - 1
509 a = src.to_a
510
511 for k in 0..size
512 if (akk = a[k][k]) == 0
513 i = k
514 begin
515 Matrix.Raise ErrNotRegular if (i += 1) > size
516 end while a[i][k] == 0
517 a[i], a[k] = a[k], a[i]
518 @rows[i], @rows[k] = @rows[k], @rows[i]
519 akk = a[k][k]
520 end
521
522 for i in 0 .. size
523 next if i == k
524 q = a[i][k] / akk
525 a[i][k] = 0
526
527 (k + 1).upto(size) do
528 |j|
529 a[i][j] -= a[k][j] * q
530 end
531 0.upto(size) do
532 |j|
533 @rows[i][j] -= @rows[k][j] * q
534 end
535 end
536
537 (k + 1).upto(size) do
538 |j|
539 a[k][j] /= akk
540 end
541 0.upto(size) do
542 |j|
543 @rows[k][j] /= akk
544 end
545 end
546 self
547 end
548 #alias reciprocal inverse
549
550 def ** (other)
551 if other.kind_of?(Integer)
552 x = self
553 if other <= 0
554 x = self.inverse
555 return Matrix.identity(self.column_size) if other == 0
556 other = -other
557 end
558 z = x
559 n = other - 1
560 while n != 0
561 while (div, mod = n.divmod(2)
562 mod == 0)
563 x = x * x
564 n = div
565 end
566 z *= x
567 n -= 1
568 end
569 z
570 elsif other.kind_of?(Float) || defined?(Rational) && other.kind_of?(Rational)
571 Matrix.Raise ErrOperationNotDefined, "**"
572 else
573 Matrix.Raise ErrOperationNotDefined, "**"
574 end
575 end
576
577 # Matrix functions
578
579 def determinant
580 return 0 unless square?
581
582 size = row_size - 1
583 a = to_a
584
585 det = 1
586 k = 0
587 begin
588 if (akk = a[k][k]) == 0
589 i = k
590 begin
591 return 0 if (i += 1) > size
592 end while a[i][k] == 0
593 a[i], a[k] = a[k], a[i]
594 akk = a[k][k]
595 det *= -1
596 end
597 (k + 1).upto(size) do
598 |i|
599 q = a[i][k] / akk
600 (k + 1).upto(size) do
601 |j|
602 a[i][j] -= a[k][j] * q
603 end
604 end
605 det *= akk
606 end while (k += 1) <= size
607 det
608 end
609 alias det determinant
610
611 def rank
612 if column_size > row_size
613 a = transpose.to_a
614 a_column_size = row_size
615 a_row_size = column_size
616 else
617 a = to_a
618 a_column_size = column_size
619 a_row_size = row_size
620 end
621 rank = 0
622 k = 0
623 begin
624 if (akk = a[k][k]) == 0
625 i = k
626 exists = true
627 begin
628 if (i += 1) > a_column_size - 1
629 exists = false
630 break
631 end
632 end while a[i][k] == 0
633 if exists
634 a[i], a[k] = a[k], a[i]
635 akk = a[k][k]
636 else
637 i = k
638 exists = true
639 begin
640 if (i += 1) > a_row_size - 1
641 exists = false
642 break
643 end
644 end while a[k][i] == 0
645 if exists
646 k.upto(a_column_size - 1) do
647 |j|
648 a[j][k], a[j][i] = a[j][i], a[j][k]
649 end
650 akk = a[k][k]
651 else
652 next
653 end
654 end
655 end
656 (k + 1).upto(a_row_size - 1) do
657 |i|
658 q = a[i][k] / akk
659 (k + 1).upto(a_column_size - 1) do
660 |j|
661 a[i][j] -= a[k][j] * q
662 end
663 end
664 rank += 1
665 end while (k += 1) <= a_column_size - 1
666 return rank
667 end
668
669 def trace
670 tr = 0
671 0.upto(column_size - 1) do
672 |i|
673 tr += @rows[i][i]
674 end
675 tr
676 end
677 alias tr trace
678
679 def transpose
680 Matrix.columns(@rows)
681 end
682 alias t transpose
683
684 # CONVERTING
685
686 def coerce(other)
687 case other
688 when Numeric
689 return Scalar.new(other), self
690 else
691 raise TypeError, "#{type} can't be coerced into #{other.type}"
692 end
693 end
694
695 def row_vectors
696 rows = (0 .. row_size - 1).collect {
697 |i|
698 row(i)
699 }
700 rows
701 end
702
703 def column_vectors
704 columns = (0 .. column_size - 1).collect {
705 |i|
706 column(i)
707 }
708 columns
709 end
710
711 def to_a
712 @rows.collect{|row| row.collect{|e| e}}
713 end
714
715 def to_f
716 collect{|e| e.to_f}
717 end
718
719 def to_i
720 collect{|e| e.to_i}
721 end
722
723 def to_r
724 collect{|e| e.to_r}
725 end
726
727 # PRINTING
728 def to_s
729 "Matrix[" + @rows.collect{
730 |row|
731 "[" + row.collect{|e| e.to_s}.join(", ") + "]"
732 }.join(", ")+"]"
733 end
734
735 def inspect
736 "Matrix"+@rows.inspect
737 end
738
739 # Private CLASS
740
741 class Scalar < Numeric
742 include ExceptionForMatrix
743
744 def initialize(value)
745 @value = value
746 end
747
748 # ARITHMETIC
749 def +(other)
750 case other
751 when Numeric
752 Scalar.new(@value + other)
753 when Vector, Matrix
754 Scalar.Raise WrongArgType, other.type, "Numeric or Scalar"
755 when Scalar
756 Scalar.new(@value + other.value)
757 else
758 x, y = other.coerce(self)
759 x + y
760 end
761 end
762
763 def -(other)
764 case other
765 when Numeric
766 Scalar.new(@value - other)
767 when Vector, Matrix
768 Scalar.Raise WrongArgType, other.type, "Numeric or Scalar"
769 when Scalar
770 Scalar.new(@value - other.value)
771 else
772 x, y = other.coerce(self)
773 x - y
774 end
775 end
776
777 def *(other)
778 case other
779 when Numeric
780 Scalar.new(@value * other)
781 when Vector, Matrix
782 other.collect{|e| @value * e}
783 else
784 x, y = other.coerce(self)
785 x * y
786 end
787 end
788
789 def / (other)
790 case other
791 when Numeric
792 Scalar.new(@value / other)
793 when Vector
794 Scalar.Raise WrongArgType, other.type, "Numeric or Scalar or Matrix"
795 when Matrix
796 self * _M.inverse
797 else
798 x, y = other.coerce(self)
799 x / y
800 end
801 end
802
803 def ** (other)
804 case other
805 when Numeric
806 Scalar.new(@value ** other)
807 when Vector
808 Scalar.Raise WrongArgType, other.type, "Numeric or Scalar or Matrix"
809 when Matrix
810 other.powered_by(self)
811 else
812 x, y = other.coerce(self)
813 x ** y
814 end
815 end
816 end
817 end
818
819 #----------------------------------------------------------------------
820 #
821 # -
822 #
823 #----------------------------------------------------------------------
824 class Vector
825 include ExceptionForMatrix
826
827 #INSTANCE CREATION
828
829 private_class_method :new
830 def Vector.[](*array)
831 new(:init_elements, array, copy = false)
832 end
833
834 def Vector.elements(array, copy = true)
835 new(:init_elements, array, copy)
836 end
837
838 def initialize(method, array, copy)
839 self.send(method, array, copy)
840 end
841
842 def init_elements(array, copy)
843 if copy
844 @elements = array.dup
845 else
846 @elements = array
847 end
848 end
849
850 # ACCSESSING
851
852 def [](i)
853 @elements[i]
854 end
855
856 def size
857 @elements.size
858 end
859
860 # ENUMRATIONS
861 def each2(v)
862 Vector.Raise ErrDimensionMismatch if size != v.size
863 0.upto(size - 1) do
864 |i|
865 yield @elements[i], v[i]
866 end
867 end
868
869 def collect2(v)
870 Vector.Raise ErrDimensionMismatch if size != v.size
871 (0 .. size - 1).collect do
872 |i|
873 yield @elements[i], v[i]
874 end
875 end
876
877 # COMPARING
878 def ==(other)
879 return false unless Vector === other
880
881 other.compare_by(@elements)
882 end
883 alias eqn? ==
884
885 def compare_by(elements)
886 @elements == elements
887 end
888
889 def clone
890 Vector.elements(@elements)
891 end
892
893 def hash
894 @elements.hash
895 end
896
897 # ARITHMETIC
898
899 def *(x) #x is matrix or number
900 case x
901 when Numeric
902 els = @elements.collect{|e| e * x}
903 Vector.elements(els, false)
904 when Matrix
905 Matrix.column_vector(self) * x
906 else
907 s, x = x.coerce(self)
908 s * x
909 end
910 end
911
912 def +(v)
913 case v
914 when Vector
915 Vector.Raise ErrDimensionMismatch if size != v.size
916 els = collect2(v) {
917 |v1, v2|
918 v1 + v2
919 }
920 Vector.elements(els, false)
921 when Matrix
922 Matrix.column_vector(self) + v
923 else
924 s, x = v.coerce(self)
925 s + x
926 end
927 end
928
929 def -(v)
930 case v
931 when Vector
932 Vector.Raise ErrDimensionMismatch if size != v.size
933 els = collect2(v) {
934 |v1, v2|
935 v1 - v2
936 }
937 Vector.elements(els, false)
938 when Matrix
939 Matrix.column_vector(self) - v
940 else
941 s, x = v.coerce(self)
942 s - x
943 end
944 end
945
946 # VECTOR FUNCTIONS
947
948 def inner_product(v)
949 Vector.Raise ErrDimensionMismatch if size != v.size
950
951 p = 0
952 each2(v) {
953 |v1, v2|
954 p += v1 * v2
955 }
956 p
957 end
958
959 def collect
960 els = @elements.collect {
961 |v|
962 yield v
963 }
964 Vector.elements(els, false)
965 end
966 alias map collect
967
968 def map2(v)
969 els = collect2(v) {
970 |v1, v2|
971 yield v1, v2
972 }
973 Vector.elements(els, false)
974 end
975
976 def r
977 v = 0
978 for e in @elements
979 v += e*e
980 end
981 return Math.sqrt(v)
982 end
983
984 # CONVERTING
985 def covector
986 Matrix.row_vector(self)
987 end
988
989 def to_a
990 @elements.dup
991 end
992
993 def to_f
994 collect{|e| e.to_f}
995 end
996
997 def to_i
998 collect{|e| e.to_i}
999 end
1000
1001 def to_r
1002 collect{|e| e.to_r}
1003 end
1004
1005 def coerce(other)
1006 case other
1007 when Numeric
1008 return Scalar.new(other), self
1009 else
1010 raise TypeError, "#{type} can't be coerced into #{other.type}"
1011 end
1012 end
1013
1014 # PRINTING
1015
1016 def to_s
1017 "Vector[" + @elements.join(", ") + "]"
1018 end
1019
1020 def inspect
1021 str = "Vector"+@elements.inspect
1022 end
1023 end