1 // ========================================================================== 2 // SeqAn - The Library for Sequence Analysis 3 // ========================================================================== 4 // Copyright (c) 2006-2015, Knut Reinert, FU Berlin 5 // All rights reserved. 6 // 7 // Redistribution and use in source and binary forms, with or without 8 // modification, are permitted provided that the following conditions are met: 9 // 10 // * Redistributions of source code must retain the above copyright 11 // notice, this list of conditions and the following disclaimer. 12 // * Redistributions in binary form must reproduce the above copyright 13 // notice, this list of conditions and the following disclaimer in the 14 // documentation and/or other materials provided with the distribution. 15 // * Neither the name of Knut Reinert or the FU Berlin nor the names of 16 // its contributors may be used to endorse or promote products derived 17 // from this software without specific prior written permission. 18 // 19 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" 20 // AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 21 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE 22 // ARE DISCLAIMED. IN NO EVENT SHALL KNUT REINERT OR THE FU BERLIN BE LIABLE 23 // FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL 24 // DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR 25 // SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER 26 // CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT 27 // LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY 28 // OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH 29 // DAMAGE. 30 // 31 // ========================================================================== 32 // Author: Andreas Gogol-Doering <andreas.doering@mdc-berlin.de> 33 // ========================================================================== 34 // Simple matrices; Used in many alignment algorithms. 35 // ========================================================================== 36 37 #ifndef SEQAN_HEADER_MATRIX_BASE_H 38 #define SEQAN_HEADER_MATRIX_BASE_H 39 40 namespace SEQAN_NAMESPACE_MAIN 41 { 42 43 ////////////////////////////////////////////////////////////////////////////// 44 45 struct NDimensional; 46 47 48 template <typename TValue, unsigned DIMENSION = 0/*typename TSpec = NDimensional*/> 49 class Matrix; 50 51 ////////////////////////////////////////////////////////////////////////////// 52 template <typename T> struct SizeArr_; 53 54 template <typename TValue, unsigned DIMENSION> 55 struct SizeArr_<Matrix<TValue, DIMENSION> > 56 { 57 typedef Matrix<TValue, DIMENSION> TMatrix_; 58 typedef typename Size<TMatrix_>::Type TSize_; 59 typedef String<TSize_> Type; 60 }; 61 62 ////////////////////////////////////////////////////////////////////////////// 63 64 template <typename TValue, unsigned DIMENSION> 65 struct Host<Matrix<TValue, DIMENSION> > 66 { 67 typedef String<TValue> Type; 68 }; 69 70 ////////////////////////////////////////////////////////////////////////////// 71 ////////////////////////////////////////////////////////////////////////////// 72 73 // TODO(holtgrew): Add more comprehensive documentation! 74 75 /*! 76 * @class Matrix 77 * @headerfile <seqan/align.h> 78 * @brief A simple n-dimensional matrix type. 79 * 80 * @signature template <typename TValue[, unsigned DIMENSION]> 81 * class Matrix; 82 * 83 * @tparam TValue Type of matrix entries. 84 * @tparam DIMENSION Dimension of the matrix. Use 0 for n-dimensional, values > 0 for a matrix with 85 * <tt>DIMENSION</tt> dimensions. Defaults to 0. 86 */ 87 88 89 template <typename TValue> 90 class Matrix<TValue, 0> 91 { 92 //____________________________________________________________________________ 93 94 public: 95 typedef typename Size<Matrix>::Type TSize; 96 typedef String<TSize> TSizeArr; 97 typedef String<TValue> THost; 98 99 TSizeArr data_lengths; //Length of every dimension 100 TSizeArr data_factors; //used for positions of dimensions in host ("size of jumps" to get to next entry of specified dimension) 101 102 Holder<THost> data_host; 103 //____________________________________________________________________________ 104 105 public: 106 Matrix() 107 { 108 create(data_host); 109 } 110 Matrix(Matrix const & other_): 111 data_lengths(other_.data_lengths), 112 data_factors(other_.data_factors), 113 data_host(other_.data_host) 114 { 115 } 116 inline Matrix const & 117 operator = (Matrix const & other_) 118 { 119 data_lengths = other_.data_lengths; 120 data_factors = other_.data_factors; 121 data_host = other_.data_host; 122 123 return *this; 124 } 125 ~Matrix() 126 { 127 } 128 //____________________________________________________________________________ 129 130 131 //____________________________________________________________________________ 132 133 inline TValue & 134 operator () (TSize x1, TSize x2) 135 { 136 return value(*this, x1, x2); 137 } 138 inline TValue & 139 operator () (TSize x1, TSize x2, TSize x3) 140 { 141 return value(*this, x1, x2, x3); 142 } 143 inline TValue & 144 operator () (TSize x1, TSize x2, TSize x3, TSize x4) 145 { 146 return value(*this, x1, x2, x3, x4); 147 } 148 149 //____________________________________________________________________________ 150 }; 151 152 153 template <typename TValue> 154 class Matrix<TValue, 2> 155 { 156 //____________________________________________________________________________ 157 158 public: 159 typedef typename Size<Matrix>::Type TSize; 160 typedef String<TSize> TSizeArr; 161 typedef String<TValue> THost; 162 163 TSizeArr data_lengths; 164 TSizeArr data_factors; 165 166 Holder<THost> data_host; 167 168 169 //____________________________________________________________________________ 170 171 public: 172 Matrix() 173 { 174 create(data_host); 175 176 //setDimension to 2 177 resize(data_lengths, 2, 0); 178 resize(data_factors, 2, 0); 179 data_factors[0] = 1; 180 } 181 Matrix(Matrix const & other_): 182 data_lengths(other_.data_lengths), 183 data_factors(other_.data_factors), 184 data_host(other_.data_host) 185 { 186 } 187 inline Matrix const & 188 operator = (Matrix const & other_) 189 { 190 data_lengths = other_.data_lengths; 191 data_factors = other_.data_factors; 192 data_host = other_.data_host; 193 194 return *this; 195 } 196 197 ~Matrix() 198 { 199 } 200 //____________________________________________________________________________ 201 202 203 //____________________________________________________________________________ 204 205 inline TValue & 206 operator () (TSize x1, TSize x2) 207 { 208 return value(*this, x1, x2); 209 } 210 211 //____________________________________________________________________________ 212 }; 213 214 template <typename TValue> 215 class Matrix<TValue, 3> 216 { 217 //____________________________________________________________________________ 218 219 public: 220 typedef typename Size<Matrix>::Type TSize; 221 typedef String<TSize> TSizeArr; 222 typedef String<TValue> THost; 223 224 TSizeArr data_lengths; 225 TSizeArr data_factors; 226 227 Holder<THost> data_host; 228 229 230 //____________________________________________________________________________ 231 232 public: 233 Matrix() 234 { 235 create(data_host); 236 237 //setDimension to 3 238 resize(data_lengths, 3, 0); 239 resize(data_factors, 3); 240 data_factors[0] = 1; 241 } 242 Matrix(Matrix const & other_): 243 data_lengths(other_.data_lengths), 244 data_factors(other_.data_factors), 245 data_host(other_.data_host) 246 { 247 } 248 inline Matrix const & 249 operator = (Matrix const & other_) 250 { 251 data_lengths = other_.data_lengths; 252 data_factors = other_.data_factors; 253 data_host = other_.data_host; 254 255 return *this; 256 } 257 258 ~Matrix() 259 { 260 } 261 //____________________________________________________________________________ 262 263 264 //____________________________________________________________________________ 265 266 inline TValue & 267 operator () (TSize x1, TSize x2, TSize x3) 268 { 269 return value(*this, x1, x2, x3); 270 } 271 272 //____________________________________________________________________________ 273 }; 274 275 template <typename TValue, unsigned DIMENSION> 276 inline typename SizeArr_<Matrix<TValue, DIMENSION> >::Type & 277 _dataLengths(Matrix<TValue, DIMENSION> & me) 278 { 279 return me.data_lengths; 280 } 281 282 template <typename TValue, unsigned DIMENSION> 283 inline typename SizeArr_<Matrix<TValue, DIMENSION> >::Type const & 284 _dataLengths(Matrix<TValue, DIMENSION> const & me) 285 { 286 return me.data_lengths; 287 } 288 289 template <typename TValue, unsigned DIMENSION> 290 inline typename SizeArr_<Matrix<TValue, DIMENSION> >::Type & 291 _dataFactors(Matrix<TValue, DIMENSION> & me) 292 { 293 return me.data_factors; 294 } 295 296 template <typename TValue, unsigned DIMENSION> 297 inline typename SizeArr_<Matrix<TValue, DIMENSION> >::Type const & 298 _dataFactors(Matrix<TValue, DIMENSION> const & me) 299 { 300 return me.data_factors; 301 } 302 303 //____________________________________________________________________________ 304 305 306 template <typename TValue, unsigned DIMENSION> 307 inline bool 308 dependent(Matrix<TValue, DIMENSION> & me) 309 { 310 return dependent(me.data_host); 311 } 312 313 //____________________________________________________________________________ 314 315 template <typename TValue, unsigned DIMENSION, typename THost> 316 inline void 317 setHost(Matrix<TValue, DIMENSION> & me, THost & host_) 318 { 319 setValue(me.data_host, host_); 320 } 321 322 //____________________________________________________________________________ 323 324 325 template <typename TValue, unsigned DIMENSION> 326 inline typename Host<Matrix<TValue, DIMENSION> >::Type & 327 host(Matrix<TValue, DIMENSION> & me) 328 { 329 return value(me.data_host); 330 } 331 332 template <typename TValue, unsigned DIMENSION> 333 inline typename Host<Matrix<TValue, DIMENSION> >::Type const & 334 host(Matrix<TValue, DIMENSION> const & me) 335 { 336 return value(me.data_host); 337 } 338 339 //____________________________________________________________________________ 340 341 template <typename TValue, unsigned DIMENSION, typename THost> 342 inline void 343 assignHost(Matrix<TValue, DIMENSION> & me, THost const & value_) 344 { 345 assignValue(me.data_host, value_); 346 } 347 348 //____________________________________________________________________________ 349 350 template <typename TValue, unsigned DIMENSION, typename THost> 351 inline void 352 moveHost(Matrix<TValue, DIMENSION> & me, THost const & value_) 353 { 354 moveValue(me.data_host, value_); 355 } 356 ////////////////////////////////////////////////////////////////////////////// 357 ////////////////////////////////////////////////////////////////////////////// 358 359 template <typename TValue, unsigned DIMENSION> 360 struct Value< Matrix<TValue, DIMENSION> > 361 { 362 typedef TValue Type; 363 }; 364 365 ////////////////////////////////////////////////////////////////////////////// 366 367 template <typename TValue, unsigned DIMENSION, typename TIteratorSpec> 368 struct Iterator< Matrix<TValue, DIMENSION>, TIteratorSpec > 369 { 370 typedef Iter<Matrix<TValue, DIMENSION>, PositionIterator> Type; 371 }; 372 373 template <typename TValue, unsigned DIMENSION, typename TIteratorSpec> 374 struct Iterator< Matrix<TValue, DIMENSION> const, TIteratorSpec > 375 { 376 typedef Iter<Matrix<TValue, DIMENSION> const, PositionIterator> Type; 377 }; 378 379 ////////////////////////////////////////////////////////////////////////////// 380 ////////////////////////////////////////////////////////////////////////////// 381 382 template <typename TValue, unsigned DIMENSION> 383 inline unsigned int 384 dimension(Matrix<TValue, DIMENSION> const & me) 385 { 386 return length(_dataLengths(me)); 387 } 388 389 ////////////////////////////////////////////////////////////////////////////// 390 391 template <typename TValue, unsigned DIMENSION> 392 inline void 393 setDimension(Matrix<TValue, DIMENSION> & me, 394 unsigned int dim_) 395 { 396 397 SEQAN_ASSERT_GT(dim_, 0u); 398 //std::cout<<"\npress enter1\n"; 399 //std::cin.get(); 400 resize(_dataLengths(me), dim_, 0); 401 402 resize(_dataFactors(me), dim_); 403 _dataFactors(me)[0] = 1; 404 } 405 406 ////////////////////////////////////////////////////////////////////////////// 407 408 template <typename TValue, unsigned DIMENSION> 409 inline typename Size<Matrix<TValue, DIMENSION> >::Type 410 length(Matrix<TValue, DIMENSION> const & me, 411 unsigned int dim_) 412 { 413 return me.data_lengths[dim_]; 414 } 415 416 template <typename TValue, unsigned DIMENSION> 417 inline typename Size<Matrix <TValue, DIMENSION> >::Type 418 length(Matrix<TValue, DIMENSION> const & me) 419 { 420 return length(host(me)); 421 } 422 423 template <typename TValue, unsigned DIMENSION> 424 inline bool empty(Matrix<TValue, DIMENSION> const & me) 425 { 426 return empty(host(me)); 427 } 428 429 ////////////////////////////////////////////////////////////////////////////// 430 431 template <typename TValue, unsigned DIMENSION, typename TLength> 432 inline void 433 setLength(Matrix<TValue, DIMENSION> & me, 434 unsigned int dim_, 435 TLength length_) 436 { 437 SEQAN_ASSERT_GT(length_, static_cast<TLength>(0)); 438 SEQAN_ASSERT_LT(dim_, dimension(me)); 439 440 typedef typename SizeArr_<Matrix<TValue, DIMENSION> >::TSize_ TSize_; 441 442 _dataLengths(me)[dim_] = static_cast<TSize_>(length_); 443 } 444 445 ////////////////////////////////////////////////////////////////////////////// 446 447 /*! 448 * @fn Matrix#resize 449 * @brief Resize the matrix and fill it with a given value or zeroes. 450 * 451 * @signature void resize(matrix[, val]); 452 * 453 * @param[in,out] matrix The Matrix to fill. 454 * @param[in] val The optional value to fill the matrix with. 455 */ 456 457 458 template <typename TValue, unsigned DIMENSION> 459 inline void 460 resize(Matrix<TValue, DIMENSION> & me) 461 { 462 typedef Matrix<TValue, DIMENSION> TMatrix; 463 typedef typename Size<TMatrix>::Type TSize; 464 465 unsigned int dimension_ = dimension(me); 466 467 SEQAN_ASSERT_GT(dimension_, 0u); 468 469 TSize factor_ = _dataFactors(me)[0] * length(me, 0); 470 for (unsigned int i = 1; (factor_ > 0) && (i < dimension_); ++i) 471 { 472 _dataFactors(me)[i] = factor_; 473 factor_ *= length(me, i); 474 } 475 476 if (factor_ > 0) 477 { 478 resize(host(me), factor_); 479 } 480 } 481 482 ////////////////////////////////////////////////////////////////////////////// 483 484 template <typename TValue, unsigned DIMENSION, typename TFillValue> 485 inline void 486 resize(Matrix<TValue, DIMENSION> & me, TFillValue myValue) //resize the matrix and fill with value 487 { 488 typedef Matrix<TValue, DIMENSION> TMatrix; 489 typedef typename Size<TMatrix>::Type TSize; 490 491 unsigned int dimension_ = dimension(me); 492 493 SEQAN_ASSERT_GT(dimension_, 0u); 494 495 TSize factor_ = _dataFactors(me)[0] * length(me, 0); 496 for (unsigned int i = 1; (factor_ > 0) && (i < dimension_); ++i) 497 { 498 _dataFactors(me)[i] = factor_; 499 factor_ *= length(me, i); 500 } 501 502 if (factor_ > 0) 503 resize(host(me), factor_, myValue); 504 } 505 506 507 ////////////////////////////////////////////////////////////////////////////// 508 509 template <typename TValue, unsigned DIMENSION, typename TPosition> 510 inline typename Position<Matrix <TValue, DIMENSION> >::Type 511 nextPosition(Matrix<TValue, DIMENSION> & me, 512 TPosition position_, 513 unsigned int dimension_) 514 { 515 return position_ + _dataFactors(me)[dimension_]; 516 } 517 518 template <typename TValue, unsigned DIMENSION, typename TPosition> 519 inline typename Position<Matrix <TValue, DIMENSION> >::Type 520 nextPosition(Matrix<TValue, DIMENSION> const & me, 521 TPosition position_, 522 unsigned int dimension_) 523 { 524 return position_ + _dataFactors(me)[dimension_]; 525 } 526 527 template <typename TValue, unsigned DIMENSION, typename TPosition> 528 inline typename Position<Matrix <TValue, DIMENSION> >::Type 529 previousPosition(Matrix<TValue, DIMENSION> & me, 530 TPosition position_, 531 unsigned int dimension_) 532 { 533 return position_ - _dataFactors(me)[dimension_]; 534 } 535 536 template <typename TValue, unsigned DIMENSION, typename TPosition> 537 inline typename Position<Matrix <TValue, DIMENSION> >::Type 538 previousPosition(Matrix<TValue, DIMENSION> const & me, 539 TPosition position_, 540 unsigned int dimension_) 541 { 542 return position_ - _dataFactors(me)[dimension_]; 543 } 544 545 ////////////////////////////////////////////////////////////////////////////// 546 547 template <typename TValue, unsigned DIMENSION, typename TPosition> 548 inline typename Size< Matrix <TValue, DIMENSION> >::Type 549 coordinate(Matrix<TValue, DIMENSION> const & me, 550 TPosition position_, 551 unsigned int dimension_) 552 { 553 SEQAN_ASSERT_LT(dimension_, dimension(me)); 554 555 if (dimension_ < dimension(me) - 1) 556 { 557 return (position_ / _dataFactors(me)[dimension_]) % _dataFactors(me)[dimension_ + 1]; 558 } 559 else 560 { 561 return position_ / _dataFactors(me)[dimension_]; 562 } 563 } 564 565 ////////////////////////////////////////////////////////////////////////////// 566 567 template <typename TValue, unsigned DIMENSION, typename TTag> 568 inline typename Iterator<Matrix <TValue, DIMENSION>, Tag<TTag> const>::Type 569 begin(Matrix<TValue, DIMENSION> & me, 570 Tag<TTag> const) 571 { 572 return typename Iterator<Matrix <TValue, DIMENSION>, Tag<TTag> const >::Type(me, 0); 573 } 574 template <typename TValue, unsigned DIMENSION, typename TTag> 575 inline typename Iterator<Matrix <TValue, DIMENSION> const, Tag<TTag> const>::Type 576 begin(Matrix<TValue, DIMENSION> const & me, 577 Tag<TTag> const) 578 { 579 return typename Iterator<Matrix <TValue, DIMENSION> const, Tag<TTag> const >::Type(me, 0); 580 } 581 582 ////////////////////////////////////////////////////////////////////////////// 583 584 template <typename TValue, unsigned DIMENSION, typename TTag> 585 inline typename Iterator<Matrix <TValue, DIMENSION>, Tag<TTag> const >::Type 586 end(Matrix<TValue, DIMENSION> & me, 587 Tag<TTag> const) 588 { 589 return typename Iterator<Matrix <TValue, DIMENSION>, Tag<TTag> const >::Type(me, length(host(me))); 590 } 591 template <typename TValue, unsigned DIMENSION, typename TTag> 592 inline typename Iterator<Matrix <TValue, DIMENSION> const, Tag<TTag> const >::Type 593 end(Matrix<TValue, DIMENSION> const & me, 594 Tag<TTag> const) 595 { 596 return typename Iterator<Matrix <TValue, DIMENSION>, Tag<TTag> const >::Type(me, length(host(me))); 597 } 598 599 ////////////////////////////////////////////////////////////////////////////// 600 601 template <typename TValue, unsigned DIMENSION, typename TPosition> 602 inline typename Reference<Matrix<TValue, DIMENSION> >::Type 603 value(Matrix<TValue, DIMENSION> & me, 604 TPosition position_) 605 { 606 return value(host(me), position_); 607 } 608 609 template <typename TValue, unsigned DIMENSION, typename TPosition> 610 inline typename Reference<Matrix<TValue, DIMENSION> const>::Type 611 value(Matrix<TValue, DIMENSION> const & me, 612 TPosition position_) 613 { 614 return value(host(me), position_); 615 } 616 617 //____________________________________________________________________________ 618 619 //two dimensional value access 620 template <typename TValue, unsigned DIMENSION, typename TOrdinate1, typename TOrdinate2> 621 inline typename Reference<Matrix<TValue, DIMENSION> >::Type 622 value(Matrix<TValue, DIMENSION> & me, 623 TOrdinate1 i1, 624 TOrdinate2 i2) 625 { 626 return value(host(me), i1 + i2 * _dataFactors(me)[1]); 627 } 628 629 template <typename TValue, unsigned DIMENSION, typename TOrdinate1, typename TOrdinate2> 630 inline typename Reference<Matrix<TValue, DIMENSION> const>::Type 631 value(Matrix<TValue, DIMENSION> const & me, 632 TOrdinate1 i1, 633 TOrdinate2 i2) 634 { 635 return value(host(me), i1 + i2 * _dataFactors(me)[1]); 636 } 637 638 //____________________________________________________________________________ 639 640 //3 dimensional value access 641 642 template <typename TValue, unsigned DIMENSION, typename TOrdinate1, typename TOrdinate2, typename TOrdinate3> 643 inline typename Reference<Matrix<TValue, DIMENSION> >::Type 644 value(Matrix<TValue, DIMENSION> & me, 645 TOrdinate1 i1, 646 TOrdinate2 i2, 647 TOrdinate3 i3) 648 { 649 return value(host(me), i1 + i2 * _dataFactors(me)[1] + i3 * _dataFactors(me)[2]); 650 } 651 652 //____________________________________________________________________________ 653 654 //4 dimensional value access 655 656 template <typename TValue, unsigned DIMENSION, typename TOrdinate1, typename TOrdinate2, typename TOrdinate3, typename TOrdinate4> 657 inline typename Reference<Matrix<TValue, DIMENSION> >::Type 658 value(Matrix<TValue, DIMENSION> & me, 659 TOrdinate1 i1, 660 TOrdinate2 i2, 661 TOrdinate3 i3, 662 TOrdinate4 i4) 663 { 664 return value(host(me), i1 + i2 * _dataFactors(me)[1] + i3 * _dataFactors(me)[2] + i4 * _dataFactors(me)[3]); 665 } 666 667 ////////////////////////////////////////////////////////////////////////////// 668 ////////////////////////////////////////////////////////////////////////////// 669 // Iterator: goNext 670 ////////////////////////////////////////////////////////////////////////////// 671 672 template <typename TValue, unsigned DIMENSION> 673 inline void 674 goNext(Iter<Matrix<TValue, DIMENSION>, PositionIterator> & me, 675 unsigned int dimension_) 676 { 677 setPosition(me, nextPosition(container(me), position(me), dimension_)); 678 } 679 680 template <typename TValue, unsigned DIMENSION> 681 inline void 682 goNext(Iter<Matrix<TValue, DIMENSION> const, PositionIterator> & me, 683 unsigned int dimension_) 684 { 685 setPosition(me, nextPosition(container(me), position(me), dimension_)); 686 } 687 688 template <typename TValue, unsigned DIMENSION> 689 inline void 690 goNext(Iter<Matrix<TValue, DIMENSION>, PositionIterator> & me) 691 { 692 goNext(me, 0); 693 } 694 695 template <typename TValue, unsigned DIMENSION> 696 inline void 697 goNext(Iter<Matrix<TValue, DIMENSION> const, PositionIterator> & me) 698 { 699 goNext(me, 0); 700 } 701 702 ////////////////////////////////////////////////////////////////////////////// 703 // Iterator: goPrevious 704 ////////////////////////////////////////////////////////////////////////////// 705 706 template <typename TValue, unsigned DIMENSION> 707 inline void 708 goPrevious(Iter< Matrix<TValue, DIMENSION>, PositionIterator > & me, 709 unsigned int dimension_) 710 { 711 setPosition(me, previousPosition(container(me), position(me), dimension_)); 712 } 713 714 template <typename TValue, unsigned DIMENSION> 715 inline void 716 goPrevious(Iter< Matrix<TValue, DIMENSION> const, PositionIterator > & me, 717 unsigned int dimension_) 718 { 719 setPosition(me, previousPosition(container(me), position(me), dimension_)); 720 } 721 722 template <typename TValue, unsigned DIMENSION> 723 inline void 724 goPrevious(Iter< Matrix<TValue, DIMENSION>, PositionIterator > & me) 725 { 726 goPrevious(me, 0); 727 } 728 729 template <typename TValue, unsigned DIMENSION> 730 inline void 731 goPrevious(Iter< Matrix<TValue, DIMENSION> const, PositionIterator > & me) 732 { 733 goPrevious(me, 0); 734 } 735 736 ////////////////////////////////////////////////////////////////////////////// 737 // goTo 738 ////////////////////////////////////////////////////////////////////////////// 739 740 template <typename TValue, unsigned DIMENSION, typename TPosition0, typename TPosition1> 741 inline void 742 goTo(Iter<Matrix<TValue, DIMENSION>, PositionIterator> & me, TPosition0 pos0, TPosition1 pos1) 743 { 744 setPosition(me, pos0 + pos1 * _dataFactors(container(me))[1]); 745 } 746 747 748 template <typename TValue, unsigned DIMENSION, typename TPosition0, typename TPosition1> 749 inline void 750 goTo(Iter<Matrix<TValue, DIMENSION> const, PositionIterator> & me, TPosition0 pos0, TPosition1 pos1) 751 { 752 setPosition(me, pos0 + pos1 * _dataFactors(container(me))[1]); 753 } 754 755 756 template <typename TValue, unsigned DIMENSION, typename TPosition0, typename TPosition1, typename TPosition2> 757 inline void 758 goTo(Iter<Matrix<TValue, DIMENSION>, PositionIterator> & me, TPosition0 pos0, TPosition1 pos1, TPosition2 pos2) 759 { 760 setPosition(me, pos0 + pos1 * _dataFactors(container(me))[1] + pos2 * _dataFactors(container(me))[2]); 761 } 762 763 764 template <typename TValue, unsigned DIMENSION, typename TPosition0, typename TPosition1, typename TPosition2> 765 inline void 766 goTo(Iter<Matrix<TValue, DIMENSION> const, PositionIterator> & me, TPosition0 pos0, TPosition1 pos1, TPosition2 pos2) 767 { 768 setPosition(me, pos0 + pos1 * _dataFactors(container(me))[1] + pos2 * _dataFactors(container(me))[2]); 769 } 770 771 ////////////////////////////////////////////////////////////////////////////// 772 // Iterator: coordinate 773 774 template <typename TValue, unsigned DIMENSION> 775 inline typename Size< Matrix<TValue, DIMENSION> >::Type 776 coordinate(Iter<Matrix<TValue, DIMENSION>, PositionIterator > & me, 777 unsigned int dimension_) 778 { 779 return coordinate(container(me), position(me), dimension_); 780 } 781 782 template <typename TValue, unsigned DIMENSION> 783 inline typename Size< Matrix<TValue, DIMENSION> >::Type 784 coordinate(Iter<Matrix<TValue, DIMENSION> const, PositionIterator > & me, 785 unsigned int dimension_) 786 { 787 return coordinate(container(me), position(me), dimension_); 788 } 789 790 /*! 791 * @fn Matrix::operator+ 792 * @brief Sum operator for the Matrix type. 793 * 794 * @signature TMatrix Matrix::operator+(lhs, rhs); 795 * 796 * @param[in] lhs First summand. 797 * @param[in] rhs Second summand. 798 * 799 * @return TMatrix The resulting matrix of same type as <tt>lhs</tt> and <tt>rhs</tt>. 800 */ 801 802 template <typename TValue,unsigned DIMENSION> 803 Matrix<TValue,DIMENSION> 804 operator + (Matrix<TValue,DIMENSION> const & matrix1,Matrix<TValue,DIMENSION> const & matrix2) 805 { 806 //the two matrices must have same dimension 807 SEQAN_ASSERT(_dataLengths(matrix1) == _dataLengths(matrix2)); 808 809 Matrix<TValue,DIMENSION> result; 810 //copy the first matrix 811 setDimension(result,length(_dataLengths(matrix1))); 812 _dataLengths(result) = _dataLengths(matrix1); 813 resize(result); 814 815 //add the matrices 816 for(unsigned int i = 0;i< length(host(result));++i) 817 { 818 value(host(result), i)=value(host(matrix1), i)+value(host(matrix2), i); 819 } 820 //Return matrix sum 821 return result; 822 } 823 824 template <typename TValue,unsigned DIMENSION> 825 Matrix<TValue,DIMENSION> 826 operator - (Matrix<TValue,DIMENSION> const & matrix1,Matrix<TValue,DIMENSION> const & matrix2) 827 { 828 //the two matrices must have same dimension 829 SEQAN_ASSERT(_dataLengths(matrix1) == _dataLengths(matrix2)); 830 831 Matrix<TValue,DIMENSION> result; 832 //resize the matrix 833 setDimension(result,length(_dataLengths(matrix1))); 834 _dataLengths(result) = _dataLengths(matrix1); 835 resize(result); 836 837 //subtract the matrices 838 for(unsigned int i = 0;i< length(host(result));++i) 839 { 840 value(host(result), i)=value(host(matrix1), i)-value(host(matrix2), i); 841 } 842 //Return matrix difference 843 return result; 844 } 845 846 template <typename TValue> 847 Matrix<TValue, 2> 848 operator * (Matrix<TValue, 2> const & matrix1, Matrix<TValue, 2> const & matrix2) 849 { 850 SEQAN_ASSERT_EQ(length(matrix1,1), length(matrix2,0)); 851 852 unsigned int nrow1=length(matrix1,0); 853 unsigned int ncol2=length(matrix2,1); 854 Matrix<TValue, 2> result; 855 //resize the matrix 856 setLength(result, 0, nrow1); 857 setLength(result, 1, ncol2); 858 resize(result,(TValue) 0); 859 860 //Matrix product 861 for(unsigned int row = 0; row < nrow1; row++) 862 { 863 for(unsigned int col = 0; col < ncol2; col++) 864 { 865 for(unsigned int colRes = 0; colRes < length(matrix1,1); colRes++) 866 { 867 value(result,row,col)+= value(host(matrix1), row + colRes * matrix1.data_factors[1])*value(host(matrix2), colRes + col * matrix2.data_factors[1]); 868 } 869 } 870 } 871 //return the matrix product 872 return result; 873 } 874 875 876 template <typename TValue> 877 Matrix<TValue, 2> 878 operator * (TValue const & scalar, Matrix<TValue, 2> const & matrix) 879 { 880 Matrix<TValue, 2> result; 881 result= matrix; 882 //scalar multiplication 883 for(unsigned int i = 0;i< length(host(result));++i) 884 { 885 value(host(result), i)*=scalar; 886 } 887 //return the matrix product 888 return result; 889 } 890 891 template <typename TValue> 892 Matrix<TValue, 2> 893 operator * (Matrix<TValue, 2> const & matrix, TValue const & scalar) 894 { 895 Matrix<TValue, 2> result; 896 result= matrix; 897 //scalar multiplication 898 for(unsigned int i = 0;i< length(host(result));++i) 899 { 900 value(host(result), i)*=scalar; 901 } 902 //return the matrix product 903 return result; 904 } 905 906 907 template <typename TValue, unsigned DIMENSION1, unsigned DIMENSION2> 908 bool 909 operator == (Matrix<TValue, DIMENSION1> const & matrix1, Matrix<TValue, DIMENSION2> const & matrix2) 910 { 911 bool result; 912 result= (matrix1.data_lengths==matrix2.data_lengths)&&(matrix1.data_factors==matrix2.data_factors)&&(value(matrix1.data_host)==value(matrix2.data_host))&&(DIMENSION1==DIMENSION2); 913 return result; 914 } 915 916 /* 917 template <typename TValue> 918 Matrix<TValue,2> 919 matricialSum(Matrix<TValue,2> &matrix1,Matrix<TValue,2> &matrix2) 920 { 921 //the two matrices must have same dimension 922 if(length(matrix1,0) != length(matrix2,0)||length(matrix1,1) != length(matrix2,1)) 923 { 924 fprintf(stderr,"Error: The two matrices have different dimensions"); 925 } 926 927 928 unsigned int nrow=length(matrix1,0); 929 unsigned int ncol=length(matrix1,1); 930 931 Matrix<TValue,2> result; 932 //resize the matrix 933 setLength(result, 0, nrow); 934 setLength(result, 1, ncol); 935 resize(result); 936 937 //add the matrices 938 for(unsigned int i = 0;i< nrow*ncol;++i) 939 { 940 value(host(result), i)=value(host(matrix1), i)+value(host(matrix2), i); 941 } 942 //Return matrix difference 943 return result; 944 945 } 946 */ 947 ////////////////////////////////////////////////////////////////////////////// 948 // _matricialDifference 949 ////////////////////////////////////////////////////////////////////////////// 950 951 /* 952 template <typename TValue> 953 inline Matrix<TValue,2> 954 matricialDifference(Matrix<TValue,2> & matrix1, Matrix<TValue,2> & matrix2) 955 { 956 //the two matrices must have same dimension 957 if(length(matrix1,0) != length(matrix2,0)||length(matrix1,1) != length(matrix2,1)) 958 { 959 fprintf(stderr,"Error: The two matrices have different dimensions"); 960 } 961 962 unsigned int nrow=length(matrix1,0); 963 unsigned int ncol=length(matrix1,1); 964 965 Matrix<TValue,2> result; 966 //resize the matrix 967 //setDimension(result, 2); 968 setLength(result, 0, nrow); 969 setLength(result, 1, ncol); 970 resize(result); 971 972 //Substract the matrices 973 for(unsigned int i1 = 0;i1< nrow;++i1) 974 { 975 for(unsigned int i2 = 0;i2<ncol;++i2) 976 { 977 value(host(result), i1 + i2 * _dataFactors(result)[1])=value(host(matrix1), i1 + i2 * _dataFactors(matrix1)[1])-value(host(matrix2), i1 + i2 * _dataFactors(matrix2)[1]); 978 } 979 980 } 981 //Return matrix difference 982 return result; 983 } 984 */ 985 986 /* 987 template <typename TValue> 988 inline Matrix<TValue, 2> 989 matricialProduct(Matrix<TValue, 2> &matrix1, 990 Matrix<TValue, 2> &matrix2) 991 { 992 //SEQAN_ASSERT_LT(dimension_, dimension(me)); 993 if(length(matrix1,1) != length(matrix2,0)) 994 { 995 fprintf(stderr,"Error: Number of columns of matrix1 is unequal to number of rows of matrix2"); 996 } 997 998 unsigned int nrow1=length(matrix1,0); 999 unsigned int ncol2=length(matrix2,1); 1000 Matrix<TValue, 2> result; 1001 //resize the matrix 1002 setLength(result, 0, nrow1); 1003 setLength(result, 1, ncol2); 1004 resize(result,(TValue) 0); 1005 1006 //Matrix product 1007 for(unsigned int row = 0; row < nrow1; row++) 1008 { 1009 for(unsigned int col = 0; col < ncol2; col++) 1010 { 1011 for(unsigned int colRes = 0; colRes < length(matrix1,1); colRes++) 1012 { 1013 value(result,row,col)+=value(matrix1, row,colRes)*value(matrix2,colRes,col); 1014 } 1015 } 1016 } 1017 //return the matrix product 1018 return result; 1019 } 1020 */ 1021 1022 /*! 1023 * @fn Matrix#transpose 1024 * @brief Tranpose a 2D Matrix. 1025 * 1026 * @signature TMatrix transpose(matrix); 1027 * 1028 * @param[in] matrix The matrix to tranpose. 1029 * @return TMatrix The resulting tranposed matrix. 1030 */ 1031 1032 template <typename TValue> 1033 Matrix<TValue,2> 1034 transpose(Matrix<TValue,2> const & matrix) 1035 { 1036 1037 unsigned int nrow=length(matrix,0); 1038 unsigned int ncol=length(matrix,1); 1039 1040 Matrix<TValue,2> result; 1041 //resize the matrix 1042 setLength(result, 0, ncol); 1043 setLength(result, 1, nrow); 1044 resize(result); 1045 1046 1047 for(unsigned int i1 = 0;i1< nrow;++i1) 1048 { 1049 for(unsigned int i2 = 0;i2<ncol;++i2) 1050 { 1051 value(host(result), i2 + i1 * _dataFactors(result)[1])=value(host(matrix), i1 + i2 * matrix.data_factors[1]); 1052 } 1053 1054 } 1055 1056 //Return transposed matrix 1057 return result; 1058 1059 } 1060 1061 1062 template < typename TValue > 1063 std::ostream& operator<<(std::ostream &out, const Matrix<TValue,2> &matrix) 1064 { 1065 for(unsigned int i1 = 0;i1< matrix.data_lengths[0];++i1) 1066 { 1067 for(unsigned int i2 = 0;i2<(matrix.data_lengths[1]-1);++i2) 1068 { 1069 out<<value(host(matrix), i1 + i2 * matrix.data_factors[1])<<"\t"; 1070 } 1071 //Last line is followd by \n instead of \t 1072 out<<value(host(matrix), i1 + (matrix.data_lengths[1]-1) *matrix.data_factors[1])<<"\n"; 1073 } 1074 1075 return out; 1076 } 1077 ////////////////////////////////////////////////////////////////////////////// 1078 /////////////////////////////////////////////////////////////// 1079 ///// READ 1080 /* 1081 * TODO(goeke) only square matrices of fixed size can be read in... 1082 */ 1083 /////////////////////////////////////////////////////////////// 1084 // template < typename TValue > 1085 // void read(FILE *file, Matrix<TValue,2> & matrix) 1086 // { 1087 // //unsigned int column_size=3; 1088 // unsigned int column_size=pow(4,5); 1089 // //read the transition matrix 1090 // setLength(matrix, 0, column_size); 1091 // setLength(matrix, 1, column_size); 1092 // resize(matrix,0.0); 1093 // for(unsigned int row=0; row<column_size; row++) 1094 // { 1095 // for(unsigned int col=0; col<column_size; col++) 1096 // { 1097 // fscanf(file,"%lf ", & value(matrix, row,col)); 1098 // } 1099 // fscanf(file,"\n"); 1100 // } 1101 // } 1102 1103 }// namespace SEQAN_NAMESPACE_MAIN 1104 1105 #endif //#ifndef SEQAN_HEADER_... 1106