PermutationMatrix.h (20748B)
1 // This file is part of Eigen, a lightweight C++ template library 2 // for linear algebra. 3 // 4 // Copyright (C) 2009 Benoit Jacob <jacob.benoit.1@gmail.com> 5 // Copyright (C) 2009-2015 Gael Guennebaud <gael.guennebaud@inria.fr> 6 // 7 // This Source Code Form is subject to the terms of the Mozilla 8 // Public License v. 2.0. If a copy of the MPL was not distributed 9 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 10 11 #ifndef EIGEN_PERMUTATIONMATRIX_H 12 #define EIGEN_PERMUTATIONMATRIX_H 13 14 namespace Eigen { 15 16 namespace internal { 17 18 enum PermPermProduct_t {PermPermProduct}; 19 20 } // end namespace internal 21 22 /** \class PermutationBase 23 * \ingroup Core_Module 24 * 25 * \brief Base class for permutations 26 * 27 * \tparam Derived the derived class 28 * 29 * This class is the base class for all expressions representing a permutation matrix, 30 * internally stored as a vector of integers. 31 * The convention followed here is that if \f$ \sigma \f$ is a permutation, the corresponding permutation matrix 32 * \f$ P_\sigma \f$ is such that if \f$ (e_1,\ldots,e_p) \f$ is the canonical basis, we have: 33 * \f[ P_\sigma(e_i) = e_{\sigma(i)}. \f] 34 * This convention ensures that for any two permutations \f$ \sigma, \tau \f$, we have: 35 * \f[ P_{\sigma\circ\tau} = P_\sigma P_\tau. \f] 36 * 37 * Permutation matrices are square and invertible. 38 * 39 * Notice that in addition to the member functions and operators listed here, there also are non-member 40 * operator* to multiply any kind of permutation object with any kind of matrix expression (MatrixBase) 41 * on either side. 42 * 43 * \sa class PermutationMatrix, class PermutationWrapper 44 */ 45 template<typename Derived> 46 class PermutationBase : public EigenBase<Derived> 47 { 48 typedef internal::traits<Derived> Traits; 49 typedef EigenBase<Derived> Base; 50 public: 51 52 #ifndef EIGEN_PARSED_BY_DOXYGEN 53 typedef typename Traits::IndicesType IndicesType; 54 enum { 55 Flags = Traits::Flags, 56 RowsAtCompileTime = Traits::RowsAtCompileTime, 57 ColsAtCompileTime = Traits::ColsAtCompileTime, 58 MaxRowsAtCompileTime = Traits::MaxRowsAtCompileTime, 59 MaxColsAtCompileTime = Traits::MaxColsAtCompileTime 60 }; 61 typedef typename Traits::StorageIndex StorageIndex; 62 typedef Matrix<StorageIndex,RowsAtCompileTime,ColsAtCompileTime,0,MaxRowsAtCompileTime,MaxColsAtCompileTime> 63 DenseMatrixType; 64 typedef PermutationMatrix<IndicesType::SizeAtCompileTime,IndicesType::MaxSizeAtCompileTime,StorageIndex> 65 PlainPermutationType; 66 typedef PlainPermutationType PlainObject; 67 using Base::derived; 68 typedef Inverse<Derived> InverseReturnType; 69 typedef void Scalar; 70 #endif 71 72 /** Copies the other permutation into *this */ 73 template<typename OtherDerived> 74 Derived& operator=(const PermutationBase<OtherDerived>& other) 75 { 76 indices() = other.indices(); 77 return derived(); 78 } 79 80 /** Assignment from the Transpositions \a tr */ 81 template<typename OtherDerived> 82 Derived& operator=(const TranspositionsBase<OtherDerived>& tr) 83 { 84 setIdentity(tr.size()); 85 for(Index k=size()-1; k>=0; --k) 86 applyTranspositionOnTheRight(k,tr.coeff(k)); 87 return derived(); 88 } 89 90 /** \returns the number of rows */ 91 inline EIGEN_DEVICE_FUNC Index rows() const { return Index(indices().size()); } 92 93 /** \returns the number of columns */ 94 inline EIGEN_DEVICE_FUNC Index cols() const { return Index(indices().size()); } 95 96 /** \returns the size of a side of the respective square matrix, i.e., the number of indices */ 97 inline EIGEN_DEVICE_FUNC Index size() const { return Index(indices().size()); } 98 99 #ifndef EIGEN_PARSED_BY_DOXYGEN 100 template<typename DenseDerived> 101 void evalTo(MatrixBase<DenseDerived>& other) const 102 { 103 other.setZero(); 104 for (Index i=0; i<rows(); ++i) 105 other.coeffRef(indices().coeff(i),i) = typename DenseDerived::Scalar(1); 106 } 107 #endif 108 109 /** \returns a Matrix object initialized from this permutation matrix. Notice that it 110 * is inefficient to return this Matrix object by value. For efficiency, favor using 111 * the Matrix constructor taking EigenBase objects. 112 */ 113 DenseMatrixType toDenseMatrix() const 114 { 115 return derived(); 116 } 117 118 /** const version of indices(). */ 119 const IndicesType& indices() const { return derived().indices(); } 120 /** \returns a reference to the stored array representing the permutation. */ 121 IndicesType& indices() { return derived().indices(); } 122 123 /** Resizes to given size. 124 */ 125 inline void resize(Index newSize) 126 { 127 indices().resize(newSize); 128 } 129 130 /** Sets *this to be the identity permutation matrix */ 131 void setIdentity() 132 { 133 StorageIndex n = StorageIndex(size()); 134 for(StorageIndex i = 0; i < n; ++i) 135 indices().coeffRef(i) = i; 136 } 137 138 /** Sets *this to be the identity permutation matrix of given size. 139 */ 140 void setIdentity(Index newSize) 141 { 142 resize(newSize); 143 setIdentity(); 144 } 145 146 /** Multiplies *this by the transposition \f$(ij)\f$ on the left. 147 * 148 * \returns a reference to *this. 149 * 150 * \warning This is much slower than applyTranspositionOnTheRight(Index,Index): 151 * this has linear complexity and requires a lot of branching. 152 * 153 * \sa applyTranspositionOnTheRight(Index,Index) 154 */ 155 Derived& applyTranspositionOnTheLeft(Index i, Index j) 156 { 157 eigen_assert(i>=0 && j>=0 && i<size() && j<size()); 158 for(Index k = 0; k < size(); ++k) 159 { 160 if(indices().coeff(k) == i) indices().coeffRef(k) = StorageIndex(j); 161 else if(indices().coeff(k) == j) indices().coeffRef(k) = StorageIndex(i); 162 } 163 return derived(); 164 } 165 166 /** Multiplies *this by the transposition \f$(ij)\f$ on the right. 167 * 168 * \returns a reference to *this. 169 * 170 * This is a fast operation, it only consists in swapping two indices. 171 * 172 * \sa applyTranspositionOnTheLeft(Index,Index) 173 */ 174 Derived& applyTranspositionOnTheRight(Index i, Index j) 175 { 176 eigen_assert(i>=0 && j>=0 && i<size() && j<size()); 177 std::swap(indices().coeffRef(i), indices().coeffRef(j)); 178 return derived(); 179 } 180 181 /** \returns the inverse permutation matrix. 182 * 183 * \note \blank \note_try_to_help_rvo 184 */ 185 inline InverseReturnType inverse() const 186 { return InverseReturnType(derived()); } 187 /** \returns the tranpose permutation matrix. 188 * 189 * \note \blank \note_try_to_help_rvo 190 */ 191 inline InverseReturnType transpose() const 192 { return InverseReturnType(derived()); } 193 194 /**** multiplication helpers to hopefully get RVO ****/ 195 196 197 #ifndef EIGEN_PARSED_BY_DOXYGEN 198 protected: 199 template<typename OtherDerived> 200 void assignTranspose(const PermutationBase<OtherDerived>& other) 201 { 202 for (Index i=0; i<rows();++i) indices().coeffRef(other.indices().coeff(i)) = i; 203 } 204 template<typename Lhs,typename Rhs> 205 void assignProduct(const Lhs& lhs, const Rhs& rhs) 206 { 207 eigen_assert(lhs.cols() == rhs.rows()); 208 for (Index i=0; i<rows();++i) indices().coeffRef(i) = lhs.indices().coeff(rhs.indices().coeff(i)); 209 } 210 #endif 211 212 public: 213 214 /** \returns the product permutation matrix. 215 * 216 * \note \blank \note_try_to_help_rvo 217 */ 218 template<typename Other> 219 inline PlainPermutationType operator*(const PermutationBase<Other>& other) const 220 { return PlainPermutationType(internal::PermPermProduct, derived(), other.derived()); } 221 222 /** \returns the product of a permutation with another inverse permutation. 223 * 224 * \note \blank \note_try_to_help_rvo 225 */ 226 template<typename Other> 227 inline PlainPermutationType operator*(const InverseImpl<Other,PermutationStorage>& other) const 228 { return PlainPermutationType(internal::PermPermProduct, *this, other.eval()); } 229 230 /** \returns the product of an inverse permutation with another permutation. 231 * 232 * \note \blank \note_try_to_help_rvo 233 */ 234 template<typename Other> friend 235 inline PlainPermutationType operator*(const InverseImpl<Other, PermutationStorage>& other, const PermutationBase& perm) 236 { return PlainPermutationType(internal::PermPermProduct, other.eval(), perm); } 237 238 /** \returns the determinant of the permutation matrix, which is either 1 or -1 depending on the parity of the permutation. 239 * 240 * This function is O(\c n) procedure allocating a buffer of \c n booleans. 241 */ 242 Index determinant() const 243 { 244 Index res = 1; 245 Index n = size(); 246 Matrix<bool,RowsAtCompileTime,1,0,MaxRowsAtCompileTime> mask(n); 247 mask.fill(false); 248 Index r = 0; 249 while(r < n) 250 { 251 // search for the next seed 252 while(r<n && mask[r]) r++; 253 if(r>=n) 254 break; 255 // we got one, let's follow it until we are back to the seed 256 Index k0 = r++; 257 mask.coeffRef(k0) = true; 258 for(Index k=indices().coeff(k0); k!=k0; k=indices().coeff(k)) 259 { 260 mask.coeffRef(k) = true; 261 res = -res; 262 } 263 } 264 return res; 265 } 266 267 protected: 268 269 }; 270 271 namespace internal { 272 template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename _StorageIndex> 273 struct traits<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime, _StorageIndex> > 274 : traits<Matrix<_StorageIndex,SizeAtCompileTime,SizeAtCompileTime,0,MaxSizeAtCompileTime,MaxSizeAtCompileTime> > 275 { 276 typedef PermutationStorage StorageKind; 277 typedef Matrix<_StorageIndex, SizeAtCompileTime, 1, 0, MaxSizeAtCompileTime, 1> IndicesType; 278 typedef _StorageIndex StorageIndex; 279 typedef void Scalar; 280 }; 281 } 282 283 /** \class PermutationMatrix 284 * \ingroup Core_Module 285 * 286 * \brief Permutation matrix 287 * 288 * \tparam SizeAtCompileTime the number of rows/cols, or Dynamic 289 * \tparam MaxSizeAtCompileTime the maximum number of rows/cols, or Dynamic. This optional parameter defaults to SizeAtCompileTime. Most of the time, you should not have to specify it. 290 * \tparam _StorageIndex the integer type of the indices 291 * 292 * This class represents a permutation matrix, internally stored as a vector of integers. 293 * 294 * \sa class PermutationBase, class PermutationWrapper, class DiagonalMatrix 295 */ 296 template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename _StorageIndex> 297 class PermutationMatrix : public PermutationBase<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime, _StorageIndex> > 298 { 299 typedef PermutationBase<PermutationMatrix> Base; 300 typedef internal::traits<PermutationMatrix> Traits; 301 public: 302 303 typedef const PermutationMatrix& Nested; 304 305 #ifndef EIGEN_PARSED_BY_DOXYGEN 306 typedef typename Traits::IndicesType IndicesType; 307 typedef typename Traits::StorageIndex StorageIndex; 308 #endif 309 310 inline PermutationMatrix() 311 {} 312 313 /** Constructs an uninitialized permutation matrix of given size. 314 */ 315 explicit inline PermutationMatrix(Index size) : m_indices(size) 316 { 317 eigen_internal_assert(size <= NumTraits<StorageIndex>::highest()); 318 } 319 320 /** Copy constructor. */ 321 template<typename OtherDerived> 322 inline PermutationMatrix(const PermutationBase<OtherDerived>& other) 323 : m_indices(other.indices()) {} 324 325 /** Generic constructor from expression of the indices. The indices 326 * array has the meaning that the permutations sends each integer i to indices[i]. 327 * 328 * \warning It is your responsibility to check that the indices array that you passes actually 329 * describes a permutation, i.e., each value between 0 and n-1 occurs exactly once, where n is the 330 * array's size. 331 */ 332 template<typename Other> 333 explicit inline PermutationMatrix(const MatrixBase<Other>& indices) : m_indices(indices) 334 {} 335 336 /** Convert the Transpositions \a tr to a permutation matrix */ 337 template<typename Other> 338 explicit PermutationMatrix(const TranspositionsBase<Other>& tr) 339 : m_indices(tr.size()) 340 { 341 *this = tr; 342 } 343 344 /** Copies the other permutation into *this */ 345 template<typename Other> 346 PermutationMatrix& operator=(const PermutationBase<Other>& other) 347 { 348 m_indices = other.indices(); 349 return *this; 350 } 351 352 /** Assignment from the Transpositions \a tr */ 353 template<typename Other> 354 PermutationMatrix& operator=(const TranspositionsBase<Other>& tr) 355 { 356 return Base::operator=(tr.derived()); 357 } 358 359 /** const version of indices(). */ 360 const IndicesType& indices() const { return m_indices; } 361 /** \returns a reference to the stored array representing the permutation. */ 362 IndicesType& indices() { return m_indices; } 363 364 365 /**** multiplication helpers to hopefully get RVO ****/ 366 367 #ifndef EIGEN_PARSED_BY_DOXYGEN 368 template<typename Other> 369 PermutationMatrix(const InverseImpl<Other,PermutationStorage>& other) 370 : m_indices(other.derived().nestedExpression().size()) 371 { 372 eigen_internal_assert(m_indices.size() <= NumTraits<StorageIndex>::highest()); 373 StorageIndex end = StorageIndex(m_indices.size()); 374 for (StorageIndex i=0; i<end;++i) 375 m_indices.coeffRef(other.derived().nestedExpression().indices().coeff(i)) = i; 376 } 377 template<typename Lhs,typename Rhs> 378 PermutationMatrix(internal::PermPermProduct_t, const Lhs& lhs, const Rhs& rhs) 379 : m_indices(lhs.indices().size()) 380 { 381 Base::assignProduct(lhs,rhs); 382 } 383 #endif 384 385 protected: 386 387 IndicesType m_indices; 388 }; 389 390 391 namespace internal { 392 template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename _StorageIndex, int _PacketAccess> 393 struct traits<Map<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime, _StorageIndex>,_PacketAccess> > 394 : traits<Matrix<_StorageIndex,SizeAtCompileTime,SizeAtCompileTime,0,MaxSizeAtCompileTime,MaxSizeAtCompileTime> > 395 { 396 typedef PermutationStorage StorageKind; 397 typedef Map<const Matrix<_StorageIndex, SizeAtCompileTime, 1, 0, MaxSizeAtCompileTime, 1>, _PacketAccess> IndicesType; 398 typedef _StorageIndex StorageIndex; 399 typedef void Scalar; 400 }; 401 } 402 403 template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename _StorageIndex, int _PacketAccess> 404 class Map<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime, _StorageIndex>,_PacketAccess> 405 : public PermutationBase<Map<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime, _StorageIndex>,_PacketAccess> > 406 { 407 typedef PermutationBase<Map> Base; 408 typedef internal::traits<Map> Traits; 409 public: 410 411 #ifndef EIGEN_PARSED_BY_DOXYGEN 412 typedef typename Traits::IndicesType IndicesType; 413 typedef typename IndicesType::Scalar StorageIndex; 414 #endif 415 416 inline Map(const StorageIndex* indicesPtr) 417 : m_indices(indicesPtr) 418 {} 419 420 inline Map(const StorageIndex* indicesPtr, Index size) 421 : m_indices(indicesPtr,size) 422 {} 423 424 /** Copies the other permutation into *this */ 425 template<typename Other> 426 Map& operator=(const PermutationBase<Other>& other) 427 { return Base::operator=(other.derived()); } 428 429 /** Assignment from the Transpositions \a tr */ 430 template<typename Other> 431 Map& operator=(const TranspositionsBase<Other>& tr) 432 { return Base::operator=(tr.derived()); } 433 434 #ifndef EIGEN_PARSED_BY_DOXYGEN 435 /** This is a special case of the templated operator=. Its purpose is to 436 * prevent a default operator= from hiding the templated operator=. 437 */ 438 Map& operator=(const Map& other) 439 { 440 m_indices = other.m_indices; 441 return *this; 442 } 443 #endif 444 445 /** const version of indices(). */ 446 const IndicesType& indices() const { return m_indices; } 447 /** \returns a reference to the stored array representing the permutation. */ 448 IndicesType& indices() { return m_indices; } 449 450 protected: 451 452 IndicesType m_indices; 453 }; 454 455 template<typename _IndicesType> class TranspositionsWrapper; 456 namespace internal { 457 template<typename _IndicesType> 458 struct traits<PermutationWrapper<_IndicesType> > 459 { 460 typedef PermutationStorage StorageKind; 461 typedef void Scalar; 462 typedef typename _IndicesType::Scalar StorageIndex; 463 typedef _IndicesType IndicesType; 464 enum { 465 RowsAtCompileTime = _IndicesType::SizeAtCompileTime, 466 ColsAtCompileTime = _IndicesType::SizeAtCompileTime, 467 MaxRowsAtCompileTime = IndicesType::MaxSizeAtCompileTime, 468 MaxColsAtCompileTime = IndicesType::MaxSizeAtCompileTime, 469 Flags = 0 470 }; 471 }; 472 } 473 474 /** \class PermutationWrapper 475 * \ingroup Core_Module 476 * 477 * \brief Class to view a vector of integers as a permutation matrix 478 * 479 * \tparam _IndicesType the type of the vector of integer (can be any compatible expression) 480 * 481 * This class allows to view any vector expression of integers as a permutation matrix. 482 * 483 * \sa class PermutationBase, class PermutationMatrix 484 */ 485 template<typename _IndicesType> 486 class PermutationWrapper : public PermutationBase<PermutationWrapper<_IndicesType> > 487 { 488 typedef PermutationBase<PermutationWrapper> Base; 489 typedef internal::traits<PermutationWrapper> Traits; 490 public: 491 492 #ifndef EIGEN_PARSED_BY_DOXYGEN 493 typedef typename Traits::IndicesType IndicesType; 494 #endif 495 496 inline PermutationWrapper(const IndicesType& indices) 497 : m_indices(indices) 498 {} 499 500 /** const version of indices(). */ 501 const typename internal::remove_all<typename IndicesType::Nested>::type& 502 indices() const { return m_indices; } 503 504 protected: 505 506 typename IndicesType::Nested m_indices; 507 }; 508 509 510 /** \returns the matrix with the permutation applied to the columns. 511 */ 512 template<typename MatrixDerived, typename PermutationDerived> 513 EIGEN_DEVICE_FUNC 514 const Product<MatrixDerived, PermutationDerived, AliasFreeProduct> 515 operator*(const MatrixBase<MatrixDerived> &matrix, 516 const PermutationBase<PermutationDerived>& permutation) 517 { 518 return Product<MatrixDerived, PermutationDerived, AliasFreeProduct> 519 (matrix.derived(), permutation.derived()); 520 } 521 522 /** \returns the matrix with the permutation applied to the rows. 523 */ 524 template<typename PermutationDerived, typename MatrixDerived> 525 EIGEN_DEVICE_FUNC 526 const Product<PermutationDerived, MatrixDerived, AliasFreeProduct> 527 operator*(const PermutationBase<PermutationDerived> &permutation, 528 const MatrixBase<MatrixDerived>& matrix) 529 { 530 return Product<PermutationDerived, MatrixDerived, AliasFreeProduct> 531 (permutation.derived(), matrix.derived()); 532 } 533 534 535 template<typename PermutationType> 536 class InverseImpl<PermutationType, PermutationStorage> 537 : public EigenBase<Inverse<PermutationType> > 538 { 539 typedef typename PermutationType::PlainPermutationType PlainPermutationType; 540 typedef internal::traits<PermutationType> PermTraits; 541 protected: 542 InverseImpl() {} 543 public: 544 typedef Inverse<PermutationType> InverseType; 545 using EigenBase<Inverse<PermutationType> >::derived; 546 547 #ifndef EIGEN_PARSED_BY_DOXYGEN 548 typedef typename PermutationType::DenseMatrixType DenseMatrixType; 549 enum { 550 RowsAtCompileTime = PermTraits::RowsAtCompileTime, 551 ColsAtCompileTime = PermTraits::ColsAtCompileTime, 552 MaxRowsAtCompileTime = PermTraits::MaxRowsAtCompileTime, 553 MaxColsAtCompileTime = PermTraits::MaxColsAtCompileTime 554 }; 555 #endif 556 557 #ifndef EIGEN_PARSED_BY_DOXYGEN 558 template<typename DenseDerived> 559 void evalTo(MatrixBase<DenseDerived>& other) const 560 { 561 other.setZero(); 562 for (Index i=0; i<derived().rows();++i) 563 other.coeffRef(i, derived().nestedExpression().indices().coeff(i)) = typename DenseDerived::Scalar(1); 564 } 565 #endif 566 567 /** \return the equivalent permutation matrix */ 568 PlainPermutationType eval() const { return derived(); } 569 570 DenseMatrixType toDenseMatrix() const { return derived(); } 571 572 /** \returns the matrix with the inverse permutation applied to the columns. 573 */ 574 template<typename OtherDerived> friend 575 const Product<OtherDerived, InverseType, AliasFreeProduct> 576 operator*(const MatrixBase<OtherDerived>& matrix, const InverseType& trPerm) 577 { 578 return Product<OtherDerived, InverseType, AliasFreeProduct>(matrix.derived(), trPerm.derived()); 579 } 580 581 /** \returns the matrix with the inverse permutation applied to the rows. 582 */ 583 template<typename OtherDerived> 584 const Product<InverseType, OtherDerived, AliasFreeProduct> 585 operator*(const MatrixBase<OtherDerived>& matrix) const 586 { 587 return Product<InverseType, OtherDerived, AliasFreeProduct>(derived(), matrix.derived()); 588 } 589 }; 590 591 template<typename Derived> 592 const PermutationWrapper<const Derived> MatrixBase<Derived>::asPermutation() const 593 { 594 return derived(); 595 } 596 597 namespace internal { 598 599 template<> struct AssignmentKind<DenseShape,PermutationShape> { typedef EigenBase2EigenBase Kind; }; 600 601 } // end namespace internal 602 603 } // end namespace Eigen 604 605 #endif // EIGEN_PERMUTATIONMATRIX_H