1 /** 2 * This module defines a matrix data structure and various operations 3 * on this data structure. All operations are @safe pure nothrow, 4 * so they can be used in any such function, and the results of any 5 * operation can be implicitly converted to an immutable type. 6 */ 7 module dstruct.matrix; 8 9 import std.traits; 10 11 // A private implementation of matrix multiplication for use in types. 12 private auto matrixMultiply(ResultType, T, U) 13 (ref ResultType result, ref const(T) left, ref const(U) right) { 14 foreach(row; 0 .. result.rowCount) { 15 foreach(column; 0 .. result.columnCount) { 16 Unqual!(typeof(result[0, 0])) value = left[row, 0] * right[0, column]; 17 18 foreach(pivot; 1 .. left.columnCount) { 19 value += left[row, pivot] * right[pivot, column]; 20 } 21 22 result[row, column] = value; 23 } 24 } 25 26 return result; 27 } 28 29 /** 30 * A matrix type. This is a 2D array of a guaranteed uniform size. 31 */ 32 struct Matrix(Number) if(isNumeric!Number) { 33 private: 34 Number[] _data; 35 size_t _rowCount; 36 size_t _columnCount; 37 public: 38 /** 39 * Create an matrix from an array of data. 40 * 41 * Params: 42 * rowCount = The number of rows for the matrix. 43 * columnCount = The number of columns for the matrix. 44 * data = The array of data for the matrix. 45 */ 46 @safe pure nothrow 47 this(size_t rowCount, size_t columnCount, immutable(Number[]) data) immutable { 48 _data = data; 49 _rowCount = rowCount; 50 _columnCount = columnCount; 51 } 52 53 // Copy-paste the constructors, because inout doesn't work with literals. 54 55 /// ditto 56 @safe pure nothrow 57 this(size_t rowCount, size_t columnCount, const(Number[]) data) const { 58 _data = data; 59 _rowCount = rowCount; 60 _columnCount = columnCount; 61 } 62 63 /// ditto 64 @safe pure nothrow 65 this(size_t rowCount, size_t columnCount, Number[] data) { 66 _data = data; 67 _rowCount = rowCount; 68 _columnCount = columnCount; 69 } 70 71 /** 72 * Create a matrix of a given size. 73 * 74 * Params: 75 * rowCount = The number of rows for the matrix. 76 * columnCount = The number of columns for the matrix. 77 */ 78 @safe pure nothrow 79 this(size_t rowCount, size_t columnCount) { 80 if (rowCount == 0 || columnCount == 0) { 81 return; 82 } 83 84 _data = new Number[](rowCount * columnCount); 85 86 _rowCount = rowCount; 87 _columnCount = columnCount; 88 } 89 90 /** 91 * Returns: A new duplicate of this matrix. 92 */ 93 @safe pure nothrow 94 Matrix!Number dup() const { 95 Matrix!Number mat; 96 97 // We can't .dup in a nothrow function, but we can do this... 98 mat._data = new Number[](_rowCount * _columnCount); 99 mat._data[] = _data[]; 100 101 mat._rowCount = _rowCount; 102 mat._columnCount = _columnCount; 103 104 return mat; 105 } 106 107 /** 108 * Returns: A new immutable duplicate of this matrix. 109 */ 110 @safe pure nothrow 111 immutable(Matrix!Number) idup() const { 112 return dup(); 113 } 114 115 /** 116 * When calling .idup on an already immutable matrix, the reference 117 * to the same immutable matrix is returned. It should be safe to 118 * share the immutable memory in this manner. 119 * 120 * Returns: A reference to this immutable matrix. 121 */ 122 @safe pure nothrow 123 immutable(Matrix!Number) idup() immutable { 124 // There's no need to copy immutable to immutable, share it! 125 return this; 126 } 127 128 unittest { 129 immutable m = immutable Matrix!int(1, 1); 130 131 // Make sure this doesn't actually duplicate. 132 assert(m._data is m._data); 133 134 auto o = Matrix!int(1, 1); 135 136 // Make sure this still does. 137 assert(o.idup._data !is o._data); 138 } 139 140 /// Returns: True if the matrix is empty. 141 @safe pure nothrow 142 @property bool empty() const { 143 return _data.length == 0; 144 } 145 146 /// Returns: The number of rows in this matrix. 147 @safe pure nothrow 148 @property size_t rowCount() const { 149 return _rowCount; 150 } 151 152 /// Returns: The number of columns in this matrix. 153 @safe pure nothrow 154 @property size_t columnCount() const { 155 return _columnCount; 156 } 157 158 /// Returns: true if the matrix is a square matrix. 159 @safe pure nothrow 160 @property bool isSquare() const { 161 return _rowCount == _columnCount; 162 } 163 164 /** 165 * Slice out a row from the matrix. Modifying this 166 * slice will modify the matrix, unless it is copied. 167 * 168 * Params: 169 * row = A row index. 170 * 171 * Returns: A slice of the row of the matrix. 172 */ 173 @trusted pure nothrow 174 inout(Number[]) opIndex(size_t row) inout 175 in { 176 assert(row <= rowCount, "row out of bounds!"); 177 } body { 178 size_t offset = row * _columnCount; 179 180 return _data[offset .. offset + _columnCount]; 181 } 182 183 /** 184 * Params: 185 * row = A row index. 186 * column = A column index. 187 * 188 * Returns: A value from the matrix 189 */ 190 @safe pure nothrow 191 ref inout(Number) opIndex(size_t row, size_t column) inout 192 in { 193 assert(column <= columnCount, "column out of bounds!"); 194 } body { 195 return _data[row * _columnCount + column]; 196 } 197 198 /** 199 * Overload for foreach(rowIndex, columnIndex, value; matrix) {} 200 */ 201 @trusted 202 int opApply(int delegate(ref size_t, ref size_t, ref Number) dg) { 203 int result = 0; 204 205 matrixLoop: foreach(row; 0 .. rowCount) { 206 size_t offset = row * columnCount; 207 208 foreach(column; 0 .. columnCount) { 209 result = dg(row, column, _data[offset + column]); 210 211 if (result) { 212 break matrixLoop; 213 } 214 } 215 } 216 217 return result; 218 } 219 220 221 /** 222 * Modify this matrix, adding/subtracting values from another matrix. 223 * 224 * Example: 225 * --- 226 * matrix += other_matrix; 227 * matrix -= yet_another_matrix; 228 * --- 229 * 230 * Params: 231 * other = Another matrix with an implicitly convertible numeric type. 232 */ 233 @safe pure nothrow 234 void opOpAssign(string op, OtherNumber) 235 (const Matrix!OtherNumber other) 236 if((op == "+" || op == "-") && is(OtherNumber : Number)) 237 in { 238 assert(this.rowCount == other.rowCount); 239 assert(this.columnCount == other.columnCount); 240 } body { 241 foreach(i; 0.._data.length) { 242 mixin(`_data[i]` ~ op ~ `= other._data[i];`); 243 } 244 } 245 246 /** 247 * Add or subtract two matrices, yielding a new matrix. 248 * 249 * Params: 250 * other = The other matrix. 251 * 252 * Returns: A new matrix. 253 */ 254 @safe pure nothrow 255 Matrix!Number opBinary(string op, OtherNumber) 256 (ref const Matrix!OtherNumber other) const 257 if((op == "+" || op == "-") && is(OtherNumber : Number)) in { 258 assert(this.rowCount == other.rowCount); 259 assert(this.columnCount == other.columnCount); 260 } out(val) { 261 assert(this.rowCount == val.rowCount); 262 assert(this.rowCount == val.rowCount); 263 } body { 264 // Copy this matrix. 265 auto result = this.dup; 266 267 mixin(`result ` ~ op ~ `= other;`); 268 269 return result; 270 } 271 272 /// ditto 273 @safe pure nothrow 274 Matrix!Number opBinary(string op, OtherNumber) 275 (const Matrix!OtherNumber other) const 276 if((op == "+" || op == "-") && is(OtherNumber : Number)) { 277 opBinary!(op, OtherNumber)(other); 278 } 279 280 /** 281 * Multiply two matrices. 282 * 283 * Given a matrix of size (m, n) and a matrix of size (o, p). 284 * This operation can only work if n == o. 285 * The resulting matrix will be size (m, p). 286 * 287 * Params: 288 * other = Another matrix 289 * 290 * Returns: The product of two matrices. 291 */ 292 @safe pure nothrow 293 Matrix!Number opBinary(string op, OtherNumber) 294 (ref const Matrix!OtherNumber other) const 295 if((op == "*") && is(OtherNumber : Number)) in { 296 assert(this.columnCount == other.rowCount); 297 } out(val) { 298 assert(val.rowCount == this.rowCount); 299 assert(val.columnCount == other.columnCount); 300 } body { 301 auto result = Matrix!Number(this.rowCount, other.columnCount); 302 303 matrixMultiply(result, this, other); 304 305 return result; 306 } 307 308 /// ditto 309 @safe pure nothrow 310 Matrix!Number opBinary(string op, OtherNumber) 311 (const Matrix!OtherNumber other) const 312 if((op == "*") && is(OtherNumber : Number)) in { 313 opBinary!(op, OtherNumber)(other); 314 } 315 316 /// ditto 317 @safe pure nothrow 318 Matrix!OtherNumber opBinary(string op, OtherNumber) 319 (ref const Matrix!OtherNumber other) const 320 if((op == "*") && !is(Number == OtherNumber) && is(Number : OtherNumber)) in { 321 assert(this.columnCount == other.rowCount); 322 } out(val) { 323 assert(val.rowCount == this.rowCount); 324 assert(val.columnCount == other.columnCount); 325 } body { 326 auto result = Matrix!OtherNumber(this.rowCount, other.columnCount); 327 328 matrixMultiply(result, this, other); 329 330 return result; 331 } 332 333 /// ditto 334 @safe pure nothrow 335 Matrix!OtherNumber opBinary(string op, OtherNumber) 336 (const Matrix!OtherNumber other) const 337 if((op == "*") && !is(Number == OtherNumber) && is(Number : OtherNumber)) { 338 opBinary!(op, OtherNumber)(other); 339 } 340 341 /** 342 * Modify this matrix with a scalar value. 343 */ 344 @safe pure nothrow 345 void opOpAssign(string op, OtherNumber)(OtherNumber other) 346 if(op != "in" && op != "~" && is(OtherNumber : Number)) { 347 foreach(i; 0.._data.length) { 348 mixin(`_data[i]` ~ op ~ `= other;`); 349 } 350 } 351 352 /** 353 * Returns: A new matrix produce by combining a matrix and a scalar value. 354 */ 355 @safe pure nothrow 356 Matrix!Number opBinary(string op, OtherNumber)(OtherNumber other) const 357 if(op != "in" && op != "~" && is(OtherNumber : Number)) { 358 // Copy this matrix. 359 auto result = this.dup; 360 361 mixin(`result ` ~ op ~ `= other;`); 362 363 return result; 364 } 365 366 /** 367 * Returns: true if two matrices are equal and have the same type. 368 */ 369 @safe pure nothrow 370 bool opEquals(ref const Matrix!Number other) const { 371 return _rowCount == other._rowCount 372 && _columnCount == other._columnCount 373 && _data == other._data; 374 } 375 376 /// ditto 377 @safe pure nothrow 378 bool opEquals(const Matrix!Number other) const { 379 return opEquals(other); 380 } 381 } 382 383 // Test basic matrix initialisation and foreach. 384 unittest { 385 size_t rowCount = 4; 386 size_t columnCount = 3; 387 int expectedValue = 42; 388 389 auto mat = Matrix!int(rowCount, columnCount, [ 390 42, 42, 42, 391 42, 42, 42, 392 42, 42, 42, 393 42, 42, 42 394 ]); 395 396 size_t expectedRow = 0; 397 size_t expectedColumn = 0; 398 399 foreach(row, column, value; mat) { 400 assert(row == expectedRow); 401 assert(column == expectedColumn); 402 assert(value == expectedValue); 403 404 if (++expectedColumn == columnCount) { 405 expectedColumn = 0; 406 ++expectedRow; 407 } 408 } 409 } 410 411 // Test matrix referencing and copying 412 unittest { 413 auto mat = Matrix!int(3, 3); 414 415 mat[0, 0] = 42; 416 417 auto normalCopy = mat.dup; 418 419 normalCopy[0, 0] = 27; 420 421 assert(mat[0, 0] == 42, "Matrix .dup created a data reference!"); 422 423 immutable immutCopy = mat.idup; 424 } 425 426 // Test modifying a matrix row externally. 427 unittest { 428 auto mat = Matrix!int(3, 3); 429 430 auto row = mat[0]; 431 432 row[0] = 3; 433 row[1] = 4; 434 row[2] = 7; 435 436 assert(mat[0] == [3, 4, 7]); 437 } 438 439 // Test immutable initialisation for a matrix 440 unittest { 441 immutable mat = immutable Matrix!int(3, 3, [ 442 1, 2, 3, 443 4, 5, 6, 444 7, 8, 9 445 ]); 446 } 447 448 // Test matrix addition/subtraction. 449 unittest { 450 void runtest(string op)() { 451 import std.stdio; 452 453 size_t rowCount = 4; 454 size_t columnCount = 4; 455 456 int leftValue = 8; 457 byte rightValue = 10; 458 459 auto left = Matrix!int(rowCount, columnCount, [ 460 8, 8, 8, 8, 461 8, 8, 8, 8, 462 8, 8, 8, 8, 463 8, 8, 8, 8, 464 ]); 465 466 auto right = Matrix!byte(rowCount, columnCount, [ 467 10, 10, 10, 10, 468 10, 10, 10, 10, 469 10, 10, 10, 10, 470 10, 10, 10, 10, 471 ]); 472 473 auto result = mixin(`left` ~ op ~ `right`); 474 auto expectedScalar = mixin(`leftValue` ~ op ~ `rightValue`); 475 476 foreach(row, column, value; result) { 477 assert(value == expectedScalar, `Matrix op failed: ` ~ op); 478 } 479 } 480 481 runtest!"+"; 482 runtest!"-"; 483 } 484 485 // Text matrix-scalar operations. 486 unittest { 487 import std.stdio; 488 489 void runtest(string op)() { 490 size_t rowCount = 2; 491 size_t columnCount = 3; 492 493 // The results for these two values are always nonzero. 494 long matrixValue = 1_234_567; 495 int scalar = 11; 496 497 auto matrix = Matrix!long(rowCount, columnCount, [ 498 matrixValue, matrixValue, matrixValue, 499 matrixValue, matrixValue, matrixValue, 500 ]); 501 502 auto result = mixin(`matrix` ~ op ~ `scalar`); 503 504 auto expectedScalar = mixin(`matrixValue` ~ op ~ `scalar`); 505 506 foreach(row, column, value; result) { 507 assert(value == expectedScalar, `Matirix scalar op failed: ` ~ op); 508 } 509 } 510 511 runtest!"+"; 512 runtest!"-"; 513 runtest!"*"; 514 runtest!"/"; 515 runtest!"%"; 516 runtest!"^^"; 517 runtest!"&"; 518 runtest!"|"; 519 runtest!"^"; 520 runtest!"<<"; 521 runtest!">>"; 522 runtest!">>>"; 523 } 524 525 unittest { 526 // Test matrix equality. 527 528 auto left = Matrix!int(3, 3, [ 529 1, 2, 3, 530 4, 5, 6, 531 7, 8, 9 532 ]); 533 534 auto right = immutable Matrix!int(3, 3, [ 535 1, 2, 3, 536 4, 5, 6, 537 7, 8, 9 538 ]); 539 540 assert(left == right); 541 } 542 543 // Test matrix multiplication 544 unittest { 545 // Let's test the Wikipedia example, why not? 546 auto left = Matrix!int(2, 3, [ 547 2, 3, 4, 548 1, 0, 0 549 ]); 550 551 auto right = Matrix!int(3, 2, [ 552 0, 1000, 553 1, 100, 554 0, 10 555 ]); 556 557 auto result = left * right; 558 559 int[][] expected = [ 560 [3, 2340], 561 [0, 1000] 562 ]; 563 564 565 foreach(i; 0..2) { 566 foreach(j; 0..2) { 567 assert(result[i, j] == expected[i][j]); 568 } 569 } 570 } 571 572 // Test matrix multiplication, with a numeric type on the 573 // left which implicitly converts to the right. 574 unittest { 575 auto left = Matrix!int(1, 1); 576 auto right = Matrix!long(1, 1); 577 578 Matrix!long result = left * right; 579 } 580 581 /** 582 * This class defines a range of rows over a matrix. 583 */ 584 struct Rows(Number) if(isNumeric!Number) { 585 private: 586 const(Number)[] _data; 587 size_t _columnCount; 588 public: 589 /** 590 * Create a new rows range for a given matrix. 591 */ 592 @safe pure nothrow 593 this(ref const Matrix!Number matrix) { 594 _data = matrix._data; 595 _columnCount = matrix._columnCount; 596 } 597 598 /// ditto 599 @safe pure nothrow 600 this(const Matrix!Number matrix) { 601 this(matrix); 602 } 603 604 /// Returns: true if the range is empty. 605 @safe pure nothrow 606 @property bool empty() const { 607 return _data.length == 0; 608 } 609 610 /// Advance to the next row. 611 @safe pure nothrow 612 void popFront() { 613 assert(!empty, "Attempted popFront on an empty Rows range!"); 614 615 _data = _data[_columnCount .. $]; 616 } 617 618 /// Returns: The current row. 619 @safe pure nothrow 620 @property const(Number[]) front() const { 621 assert(!empty, "Cannot get the front of an empty Rows range!"); 622 623 return this[0]; 624 } 625 626 /// Save a copy of this range. 627 @safe pure nothrow 628 Rows!Number save() const { 629 return this; 630 } 631 632 /// Retreat a row backwards. 633 @safe pure nothrow 634 void popBack() { 635 assert(!empty, "Attempted popBack on an empty Rows range!"); 636 637 _data = _data[0 .. $ - _columnCount]; 638 } 639 640 /// Returns: The row at the end of the range. 641 @safe pure nothrow 642 @property const(Number[]) back() const { 643 assert(!empty, "Cannot get the back of an empty Rows range!"); 644 645 return this[$ - 1]; 646 } 647 648 /** 649 * Params: 650 * index = An index for a row in the range. 651 * 652 * Returns: A row at an index in the range. 653 */ 654 @safe pure nothrow 655 @property const(Number[]) opIndex(size_t index) const in { 656 assert(index >= 0, "Negative index given to Rows opIndex!"); 657 assert(index < length, "Out of bounds index given to Rows opIndex!"); 658 } body { 659 size_t offset = index * _columnCount; 660 661 return _data[offset .. offset + _columnCount]; 662 } 663 664 /// Returns: The current length of the range. 665 @safe pure nothrow 666 @property size_t length() const { 667 if (_data.length == 0) { 668 return 0; 669 } 670 671 return _data.length / _columnCount; 672 } 673 674 /// ditto 675 @safe pure nothrow 676 @property size_t opDollar() const { 677 return length; 678 } 679 } 680 681 /** 682 * Returns: A range through a matrix's rows. 683 */ 684 @safe pure nothrow 685 Rows!Number rows(Number)(Matrix!Number matrix) { 686 return typeof(return)(matrix); 687 } 688 689 unittest { 690 auto mat = Matrix!int(3, 3, [ 691 1, 2, 3, 692 4, 5, 6, 693 7, 8, 9 694 ]); 695 696 auto expected = [ 697 [1, 2, 3], 698 [4, 5, 6], 699 [7, 8, 9] 700 ]; 701 702 size_t rowIndex = 0; 703 704 // Test InputRange stuff 705 for (auto range = mat.rows; !range.empty; range.popFront) { 706 assert(range.front == expected[rowIndex++]); 707 } 708 709 // Test ForwardRange 710 711 auto range1 = mat.rows; 712 713 range1.popFront; 714 range1.popBack; 715 716 auto range2 = range1.save; 717 718 range1.popFront; 719 720 assert(range2.front == [4, 5, 6]); 721 722 rowIndex = 3; 723 724 // Test BidirectionalRange 725 for (auto range = mat.rows; !range.empty; range.popBack) { 726 assert(range.back == expected[--rowIndex]); 727 } 728 729 // Test RandomAccessRange 730 731 auto range3 = mat.rows; 732 733 range3.popFront; 734 range3.popBack; 735 736 assert(range3.length == 1); 737 assert(range3[0] == [4, 5, 6]); 738 } 739 740 // Test 0 size Matrix rows 741 unittest { 742 Matrix!int mat; 743 744 assert(mat.rows.length == 0); 745 } 746 747 /** 748 * A static matrix type. This is a value matrix value type created directly 749 * on the stack. 750 */ 751 struct Matrix(Number, size_t _rowCount, size_t _columnCount) 752 if(isNumeric!Number && _rowCount > 0 && _columnCount > 0) { 753 /// The number of rows in this matrix. 754 enum rowCount = _rowCount; 755 /// The number of columns in this matrix. 756 enum columnCount = _columnCount; 757 /// true if this matrix is a zero-sized matrix. 758 enum empty = false; 759 /// true if this matrix is a square matrix. 760 enum isSquare = rowCount == columnCount; 761 762 /// The data backing this matrix. 763 Number[columnCount][rowCount] array2D; 764 765 alias array2D this; 766 767 /** 768 * Construct this matrix from a 2 dimensional static array. 769 * 770 * Params: 771 * array2D = A 2 dimension array of the same size. 772 */ 773 @safe pure nothrow 774 this(ref const(Number[columnCount][rowCount]) data) inout { 775 array2D = data; 776 } 777 778 /// ditto 779 @safe pure nothrow 780 this(const(Number[columnCount][rowCount]) data) inout { 781 array2D = data; 782 } 783 784 /** 785 * Construct this matrix directly from a series of numbers. 786 * This constructor is designed to be executed at compile time. 787 * 788 * Params: 789 * numbers... = A series of numbers to initialise the matrix with. 790 */ 791 @safe pure nothrow 792 this(Number[rowCount * columnCount] numbers...) { 793 foreach(row; 0 .. rowCount) { 794 foreach(column; 0 .. columnCount) { 795 array2D[row][column] = numbers[row * columnCount + column]; 796 } 797 } 798 } 799 800 /** 801 * Returns: A reference to this matrix's data as a 1D array. 802 */ 803 @trusted pure nothrow 804 @property 805 ref inout(Number[rowCount * columnCount]) array1D() inout { 806 return (cast(Number*)array2D.ptr)[0 .. rowCount * columnCount]; 807 } 808 809 // Even with alias this, we still need this overload. 810 /** 811 * Params: 812 * row = A row index. 813 * 814 * Returns: A row from the matrix. 815 */ 816 @safe pure nothrow 817 ref inout(Number[columnCount]) opIndex(size_t row) inout { 818 return array2D[row]; 819 } 820 821 /** 822 * Params: 823 * row = A row index. 824 * column = A column index. 825 * 826 * Returns: A value from the matrix 827 */ 828 @safe pure nothrow 829 ref inout(Number) opIndex(size_t row, size_t column) inout { 830 return array2D[row][column]; 831 } 832 833 /** 834 * Overload for foreach(rowIndex, columnIndex, value; matrix) {} 835 */ 836 @trusted 837 int opApply(int delegate(ref size_t, ref size_t, ref Number) dg) { 838 int result = 0; 839 840 matrixLoop: foreach(row, rowArray; array2D) { 841 foreach(column, value; rowArray) { 842 result = dg(row, column, value); 843 844 if (result) { 845 break matrixLoop; 846 } 847 } 848 } 849 850 return result; 851 } 852 853 /** 854 * Modify this matrix, adding/subtracting values from another matrix. 855 * 856 * Example: 857 * --- 858 * matrix += other_matrix; 859 * matrix -= yet_another_matrix; 860 * --- 861 * 862 * Params: 863 * other = Another matrix with an implicitly convertible numeric type. 864 */ 865 @safe pure nothrow 866 void opOpAssign(string op, OtherNumber) 867 (ref const(Matrix!(OtherNumber, rowCount, columnCount)) other) 868 if((op == "+" || op == "-") && is(OtherNumber : Number)) { 869 foreach(i; 0 .. rowCount) { 870 foreach(j; 0 .. columnCount) { 871 mixin(`array2D[i][j]` ~ op ~ `= other.array2D[i][j];`); 872 } 873 } 874 } 875 876 /// ditto 877 @safe pure nothrow 878 void opOpAssign(string op, OtherNumber) 879 (const(Matrix!(OtherNumber, rowCount, columnCount)) other) 880 if((op == "+" || op == "-") && is(OtherNumber : Number)) { 881 opOpAssign(other); 882 } 883 884 /** 885 * Add or subtract two matrices, yielding a new matrix. 886 * 887 * Params: 888 * other = The other matrix. 889 * 890 * Returns: A new matrix. 891 */ 892 @safe pure nothrow 893 Matrix!(Number, rowCount, columnCount) opBinary(string op, OtherNumber) 894 (ref const(Matrix!(OtherNumber, rowCount, columnCount)) other) 895 if((op == "+" || op == "-") && is(OtherNumber : Number)) { 896 // Copy this matrix. 897 typeof(return) result = this; 898 899 mixin(`result ` ~ op ~ `= other;`); 900 901 return result; 902 } 903 904 /// ditto 905 @safe pure nothrow 906 Matrix!(Number, rowCount, columnCount) opBinary(string op, OtherNumber) 907 (const(Matrix!(OtherNumber, rowCount, columnCount)) other) 908 if((op == "+" || op == "-") && is(OtherNumber : Number)) { 909 return opBinary(other); 910 } 911 912 /** 913 * Multiply two matrices. 914 * 915 * Given a matrix of size (m, n) and a matrix of size (o, p). 916 * This operation can only work if n == o. 917 * The resulting matrix will be size (m, p). 918 * 919 * Params: 920 * other = Another matrix 921 * 922 * Returns: The product of two matrices. 923 */ 924 @safe pure nothrow 925 Matrix!(Number, rowCount, otherColumnCount) 926 opBinary(string op, OtherNumber, size_t otherRowCount, size_t otherColumnCount) 927 (ref const(Matrix!(OtherNumber, otherRowCount, otherColumnCount)) other) const 928 if((op == "*") && (columnCount == otherRowCount) && is(OtherNumber : Number)) { 929 typeof(return) result; 930 931 matrixMultiply(result, this, other); 932 933 return result; 934 } 935 936 /// ditto 937 @safe pure nothrow 938 Matrix!(Number, rowCount, otherColumnCount) 939 opBinary(string op, OtherNumber, size_t otherRowCount, size_t otherColumnCount) 940 (const(Matrix!(OtherNumber, otherRowCount, otherColumnCount)) other) const 941 if((op == "*") && (columnCount == otherRowCount) && is(OtherNumber : Number)) in { 942 return opBinary!(op, OtherNumber, otherRowCount, otherColumnCount)(other); 943 } 944 945 @safe pure nothrow 946 Matrix!(OtherNumber, rowCount, otherColumnCount) 947 opBinary(string op, OtherNumber, size_t otherRowCount, size_t otherColumnCount) 948 (ref const(Matrix!(OtherNumber, otherRowCount, otherColumnCount)) other) const 949 if((op == "*") && (columnCount == otherRowCount) && !is(Number == OtherNumber) && is(Number : OtherNumber)) { 950 typeof(return) result; 951 952 matrixMultiply(result, this, other); 953 954 return result; 955 } 956 957 @safe pure nothrow 958 Matrix!(OtherNumber, rowCount, otherColumnCount) 959 opBinary(string op, OtherNumber, size_t otherRowCount, size_t otherColumnCount) 960 (const(Matrix!(OtherNumber, otherRowCount, otherColumnCount)) other) const 961 if((op == "*") && (columnCount == otherRowCount) && !is(Number == OtherNumber) && is(Number : OtherNumber)) { 962 return opBinary!(op, OtherNumber, otherRowCount, otherColumnCount)(other); 963 } 964 } 965 966 // 0 size matrices are special case. 967 968 /// ditto 969 struct Matrix(Number, size_t _rowCount, size_t _columnCount) 970 if(isNumeric!Number && _rowCount == 0 && _columnCount == 0) { 971 /// The number of rows in this matrix. 972 enum rowCount = _rowCount; 973 /// The number of columns in this matrix. 974 enum columnCount = _columnCount; 975 /// True if this matrix is a zero-sized matrix. 976 enum empty = true; 977 /// true if this matrix is a square matrix. 978 enum isSquare = true; 979 } 980 981 /** 982 * Alias all (M, 0), (0, N) size static matrices into one single zero-sized 983 * type of size (0, 0). 984 */ 985 template Matrix(Number, size_t rowCount, size_t columnCount) 986 if ((rowCount > 0 && columnCount == 0) || (rowCount == 0 && columnCount > 0)) { 987 alias Matrix = Matrix!(Number, 0, 0); 988 } 989 990 // Test copy constructor for 2D arrays. 991 unittest { 992 int[3][2] data = [ 993 [1, 2, 3], 994 [4, 5, 6], 995 ]; 996 997 Matrix!(int, 2, 3) matrix = data; 998 999 assert(data == matrix); 1000 } 1001 1002 // Test move constructor for 2D arrays. 1003 unittest { 1004 auto matrix = Matrix!(int, 3, 3)([ 1005 [1, 2, 3], 1006 [4, 5, 6], 1007 [7, 8, 9] 1008 ]); 1009 } 1010 1011 // Test immutable too. 1012 unittest { 1013 immutable(int[3][3]) data = [ 1014 [1, 2, 3], 1015 [4, 5, 6], 1016 [7, 8, 9] 1017 ]; 1018 1019 Matrix!(int, 3, 3) matrix = data; 1020 1021 assert(data == matrix); 1022 } 1023 1024 unittest { 1025 int[3][3] data = [ 1026 [1, 2, 3], 1027 [4, 5, 6], 1028 [7, 8, 9] 1029 ]; 1030 1031 immutable(Matrix!(int, 3, 3,)) matrix = data; 1032 1033 assert(data == matrix); 1034 } 1035 1036 // Test 1D array matrix slicing. 1037 unittest { 1038 auto matrix = Matrix!(int, 3, 3)([ 1039 [1, 2, 3], 1040 [4, 5, 6], 1041 [7, 8, 9] 1042 ]); 1043 1044 assert(matrix.array1D == [1, 2, 3, 4, 5, 6, 7, 8, 9]); 1045 1046 matrix.array1D[0] = 347; 1047 1048 assert(matrix[0][0] == 347); 1049 } 1050 1051 // Test that copy semantics work property for 1D arrays. 1052 unittest { 1053 auto matrix = Matrix!(int, 3, 3)([ 1054 [1, 2, 3], 1055 [4, 5, 6], 1056 [7, 8, 9] 1057 ]); 1058 1059 auto arr = matrix.array1D; 1060 1061 arr[0] = 347; 1062 1063 // The value should not have changed. 1064 assert(matrix[0][0] == 1); 1065 } 1066 1067 // Test that reference semantics work property for 1D arrays. 1068 unittest { 1069 auto matrix = Matrix!(int, 3, 3)([ 1070 [1, 2, 3], 1071 [4, 5, 6], 1072 [7, 8, 9] 1073 ]); 1074 1075 void foo(ref int[9] arr) { 1076 arr[0] = 347; 1077 } 1078 1079 foo(matrix.array1D); 1080 1081 // The should have changed. 1082 assert(matrix[0][0] == 347); 1083 } 1084 1085 // Test const and immutable 1D array, just in case. 1086 unittest { 1087 const matrix = Matrix!(int, 3, 3)([ 1088 [1, 2, 3], 1089 [4, 5, 6], 1090 [7, 8, 9] 1091 ]); 1092 1093 assert(matrix.array1D == [1, 2, 3, 4, 5, 6, 7, 8, 9]); 1094 assert(is(typeof(matrix.array1D) == const int[9])); 1095 } 1096 1097 unittest { 1098 immutable matrix = Matrix!(int, 3, 3)([ 1099 [1, 2, 3], 1100 [4, 5, 6], 1101 [7, 8, 9] 1102 ]); 1103 1104 assert(matrix.array1D == [1, 2, 3, 4, 5, 6, 7, 8, 9]); 1105 assert(is(typeof(matrix.array1D) == immutable int[9])); 1106 } 1107 1108 // Test compile time init with numbers. 1109 unittest { 1110 enum matrix = Matrix!(int, 3, 3)( 1111 1, 2, 3, 1112 4, 5, 6, 1113 7, 8, 9 1114 ); 1115 1116 assert(matrix.array2D[0][0] == 1); 1117 assert(matrix.array2D[0][1] == 2); 1118 assert(matrix.array2D[0][2] == 3); 1119 assert(matrix.array2D[1][0] == 4); 1120 assert(matrix.array2D[1][1] == 5); 1121 assert(matrix.array2D[1][2] == 6); 1122 assert(matrix.array2D[2][0] == 7); 1123 assert(matrix.array2D[2][1] == 8); 1124 assert(matrix.array2D[2][2] == 9); 1125 } 1126 1127 // Test zero sized matrices 1128 unittest { 1129 Matrix!(int, 0, 0) nothing; 1130 Matrix!(int, 1, 1) scalar; 1131 1132 assert(nothing.empty); 1133 assert(!scalar.empty); 1134 1135 Matrix!(int, 0, 1) noRows; 1136 Matrix!(int, 1, 0) noColumns; 1137 1138 assert(is(typeof(nothing) == typeof(noRows))); 1139 assert(is(typeof(noRows) == typeof(noColumns))); 1140 } 1141 1142 // Test index and assignment for cells. 1143 unittest { 1144 Matrix!(int, 3, 3) matrix; 1145 1146 assert(matrix[0, 0] == 0); 1147 1148 matrix[0, 0] = 3; 1149 1150 assert(matrix[0, 0] == 3); 1151 1152 matrix[0, 0] *= 3; 1153 1154 assert(matrix[0, 0] == 9); 1155 } 1156 1157 /** 1158 * Transpose (flip) a matrix. 1159 * 1160 * Params: 1161 * matrix = The matrix to produce a transpose for. 1162 * 1163 * Returns: A new matrix which is the transpose of the given matrix. 1164 */ 1165 @safe pure nothrow 1166 Matrix!Number transpose(Number)(const Matrix!Number matrix) 1167 out(val) { 1168 assert(matrix.columnCount == val.rowCount); 1169 assert(matrix.rowCount == val.columnCount); 1170 } body { 1171 auto result = typeof(return)(matrix.columnCount, matrix.rowCount); 1172 1173 foreach(row; 0 .. matrix.rowCount) { 1174 foreach(col; 0 .. matrix.columnCount) { 1175 result[col, row] = matrix[row, col]; 1176 } 1177 } 1178 1179 return result; 1180 } 1181 1182 /// ditto 1183 @safe pure nothrow 1184 Matrix!(Number, columnCount, rowCount) 1185 transpose(Number, size_t rowCount, size_t columnCount) 1186 (ref const Matrix!(Number, rowCount, columnCount) matrix) { 1187 typeof(return) result; 1188 1189 foreach(row; 0 .. matrix.rowCount) { 1190 foreach(col; 0 .. matrix.columnCount) { 1191 result[col, row] = matrix[row, col]; 1192 } 1193 } 1194 1195 return result; 1196 } 1197 1198 /// ditto 1199 @safe pure nothrow 1200 Matrix!(Number, columnCount, rowCount) 1201 transpose(Number, size_t rowCount, size_t columnCount) 1202 (const Matrix!(Number, rowCount, columnCount) matrix) { 1203 return transpose(matrix); 1204 } 1205 1206 unittest { 1207 auto matrix = Matrix!int(2, 3, [ 1208 1, 2, 3, 1209 0, -6, 7 1210 ]); 1211 1212 int[] expected = [1, 0, 2, -6, 3, 7]; 1213 1214 auto result = matrix.transpose; 1215 1216 assert(result.rowCount == matrix.columnCount); 1217 assert(result.columnCount == matrix.rowCount); 1218 assert(result._data == expected); 1219 } 1220 1221 unittest { 1222 // When transposed twice, we should get the same matrix. 1223 auto matrix = Matrix!int(2, 3, [ 1224 1, 2, 3, 1225 0, -6, 7 1226 ]); 1227 1228 assert(matrix == matrix.transpose.transpose); 1229 } 1230 1231 unittest { 1232 auto matrix = Matrix!(int, 2, 3)( 1233 1, 2, 3, 1234 0, -6, 7 1235 ); 1236 1237 assert(matrix == matrix.transpose.transpose); 1238 } 1239 1240 // Test foreach on static matrices. 1241 unittest { 1242 auto matrix = Matrix!(int, 2, 3)( 1243 1, 2, 3, 1244 0, -6, 7 1245 ); 1246 1247 foreach(row, col, value; matrix) { 1248 if (row == 0) { 1249 if (col == 0) { 1250 assert(value == 1); 1251 } else if (col == 1) { 1252 assert(value == 2); 1253 } else { 1254 assert(value == 3); 1255 } 1256 } else { 1257 if (col == 0) { 1258 assert(value == 0); 1259 } else if (col == 1) { 1260 assert(value == -6); 1261 } else { 1262 assert(value == 7); 1263 } 1264 } 1265 } 1266 } 1267 1268 // Test binary modifying operations on static matrices. 1269 unittest { 1270 auto left = Matrix!(int, 2, 3)( 1271 1, 2, 3, 1272 4, 5, 6 1273 ); 1274 1275 auto right = Matrix!(int, 2, 3)( 1276 1, 2, 3, 1277 4, 5, 6 1278 ); 1279 1280 left -= right; 1281 1282 foreach(row, col, value; left) { 1283 assert(value == 0); 1284 } 1285 1286 left -= right; 1287 left += right; 1288 1289 foreach(row, col, value; left) { 1290 assert(value == 0); 1291 } 1292 } 1293 1294 // Test binary copying operations on static matrices. 1295 unittest { 1296 auto left = Matrix!(int, 2, 3)( 1297 1, 2, 3, 1298 4, 5, 6 1299 ); 1300 1301 auto right = Matrix!(int, 2, 3)( 1302 1, 2, 3, 1303 4, 5, 6 1304 ); 1305 1306 auto newMatrix = left - right; 1307 1308 foreach(row, col, value; newMatrix) { 1309 assert(value == 0); 1310 } 1311 1312 auto finalMatrix = newMatrix + left; 1313 1314 finalMatrix -= left; 1315 1316 foreach(row, col, value; finalMatrix) { 1317 assert(value == 0); 1318 } 1319 } 1320 1321 // Test matrix multiplication on static matrices 1322 unittest { 1323 auto left = Matrix!(int, 2, 3)( 1324 2, 3, 4, 1325 1, 0, 0 1326 ); 1327 1328 auto right = Matrix!(int, 3, 2)( 1329 0, 1000, 1330 1, 100, 1331 0, 10 1332 ); 1333 1334 auto result = left * right; 1335 1336 int[][] expected = [ 1337 [3, 2340], 1338 [0, 1000] 1339 ]; 1340 1341 foreach(i; 0..2) { 1342 foreach(j; 0..2) { 1343 assert(result[i, j] == expected[i][j]); 1344 } 1345 } 1346 } 1347 1348 // Test matrix multiplication, with a numeric type on the 1349 // left which implicitly converts to the right. 1350 unittest { 1351 Matrix!(int, 1, 1) left; 1352 Matrix!(long, 1, 1) right; 1353 1354 Matrix!(long, 1, 1) result = left * right; 1355 }