PermutationMatrix.h
Go to the documentation of this file.
00001 // This file is part of Eigen, a lightweight C++ template library
00002 // for linear algebra.
00003 //
00004 // Copyright (C) 2009 Benoit Jacob <jacob.benoit.1@gmail.com>
00005 // Copyright (C) 2009-2011 Gael Guennebaud <gael.guennebaud@inria.fr>
00006 //
00007 // Eigen is free software; you can redistribute it and/or
00008 // modify it under the terms of the GNU Lesser General Public
00009 // License as published by the Free Software Foundation; either
00010 // version 3 of the License, or (at your option) any later version.
00011 //
00012 // Alternatively, you can redistribute it and/or
00013 // modify it under the terms of the GNU General Public License as
00014 // published by the Free Software Foundation; either version 2 of
00015 // the License, or (at your option) any later version.
00016 //
00017 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
00018 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
00019 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
00020 // GNU General Public License for more details.
00021 //
00022 // You should have received a copy of the GNU Lesser General Public
00023 // License and a copy of the GNU General Public License along with
00024 // Eigen. If not, see <http://www.gnu.org/licenses/>.
00025 
00026 #ifndef EIGEN_PERMUTATIONMATRIX_H
00027 #define EIGEN_PERMUTATIONMATRIX_H
00028 
00029 namespace Eigen { 
00030 
00031 template<int RowCol,typename IndicesType,typename MatrixType, typename StorageKind> class PermutedImpl;
00032 
00057 namespace internal {
00058 
00059 template<typename PermutationType, typename MatrixType, int Side, bool Transposed=false>
00060 struct permut_matrix_product_retval;
00061 template<typename PermutationType, typename MatrixType, int Side, bool Transposed=false>
00062 struct permut_sparsematrix_product_retval;
00063 enum PermPermProduct_t {PermPermProduct};
00064 
00065 } // end namespace internal
00066 
00067 template<typename Derived>
00068 class PermutationBase : public EigenBase<Derived>
00069 {
00070     typedef internal::traits<Derived> Traits;
00071     typedef EigenBase<Derived> Base;
00072   public:
00073 
00074     #ifndef EIGEN_PARSED_BY_DOXYGEN
00075     typedef typename Traits::IndicesType IndicesType;
00076     enum {
00077       Flags = Traits::Flags,
00078       CoeffReadCost = Traits::CoeffReadCost,
00079       RowsAtCompileTime = Traits::RowsAtCompileTime,
00080       ColsAtCompileTime = Traits::ColsAtCompileTime,
00081       MaxRowsAtCompileTime = Traits::MaxRowsAtCompileTime,
00082       MaxColsAtCompileTime = Traits::MaxColsAtCompileTime
00083     };
00084     typedef typename Traits::Scalar Scalar;
00085     typedef typename Traits::Index Index;
00086     typedef Matrix<Scalar,RowsAtCompileTime,ColsAtCompileTime,0,MaxRowsAtCompileTime,MaxColsAtCompileTime>
00087             DenseMatrixType;
00088     typedef PermutationMatrix<IndicesType::SizeAtCompileTime,IndicesType::MaxSizeAtCompileTime,Index>
00089             PlainPermutationType;
00090     using Base::derived;
00091     #endif
00092 
00094     template<typename OtherDerived>
00095     Derived& operator=(const PermutationBase<OtherDerived>& other)
00096     {
00097       indices() = other.indices();
00098       return derived();
00099     }
00100 
00102     template<typename OtherDerived>
00103     Derived& operator=(const TranspositionsBase<OtherDerived>& tr)
00104     {
00105       setIdentity(tr.size());
00106       for(Index k=size()-1; k>=0; --k)
00107         applyTranspositionOnTheRight(k,tr.coeff(k));
00108       return derived();
00109     }
00110 
00111     #ifndef EIGEN_PARSED_BY_DOXYGEN
00112 
00115     Derived& operator=(const PermutationBase& other)
00116     {
00117       indices() = other.indices();
00118       return derived();
00119     }
00120     #endif
00121 
00123     inline Index rows() const { return indices().size(); }
00124 
00126     inline Index cols() const { return indices().size(); }
00127 
00129     inline Index size() const { return indices().size(); }
00130 
00131     #ifndef EIGEN_PARSED_BY_DOXYGEN
00132     template<typename DenseDerived>
00133     void evalTo(MatrixBase<DenseDerived>& other) const
00134     {
00135       other.setZero();
00136       for (int i=0; i<rows();++i)
00137         other.coeffRef(indices().coeff(i),i) = typename DenseDerived::Scalar(1);
00138     }
00139     #endif
00140 
00145     DenseMatrixType toDenseMatrix() const
00146     {
00147       return derived();
00148     }
00149 
00151     const IndicesType& indices() const { return derived().indices(); }
00153     IndicesType& indices() { return derived().indices(); }
00154 
00157     inline void resize(Index size)
00158     {
00159       indices().resize(size);
00160     }
00161 
00163     void setIdentity()
00164     {
00165       for(Index i = 0; i < size(); ++i)
00166         indices().coeffRef(i) = i;
00167     }
00168 
00171     void setIdentity(Index size)
00172     {
00173       resize(size);
00174       setIdentity();
00175     }
00176 
00186     Derived& applyTranspositionOnTheLeft(Index i, Index j)
00187     {
00188       eigen_assert(i>=0 && j>=0 && i<size() && j<size());
00189       for(Index k = 0; k < size(); ++k)
00190       {
00191         if(indices().coeff(k) == i) indices().coeffRef(k) = j;
00192         else if(indices().coeff(k) == j) indices().coeffRef(k) = i;
00193       }
00194       return derived();
00195     }
00196 
00205     Derived& applyTranspositionOnTheRight(Index i, Index j)
00206     {
00207       eigen_assert(i>=0 && j>=0 && i<size() && j<size());
00208       std::swap(indices().coeffRef(i), indices().coeffRef(j));
00209       return derived();
00210     }
00211 
00216     inline Transpose<PermutationBase> inverse() const
00217     { return derived(); }
00222     inline Transpose<PermutationBase> transpose() const
00223     { return derived(); }
00224 
00225     /**** multiplication helpers to hopefully get RVO ****/
00226 
00227   
00228 #ifndef EIGEN_PARSED_BY_DOXYGEN
00229   protected:
00230     template<typename OtherDerived>
00231     void assignTranspose(const PermutationBase<OtherDerived>& other)
00232     {
00233       for (int i=0; i<rows();++i) indices().coeffRef(other.indices().coeff(i)) = i;
00234     }
00235     template<typename Lhs,typename Rhs>
00236     void assignProduct(const Lhs& lhs, const Rhs& rhs)
00237     {
00238       eigen_assert(lhs.cols() == rhs.rows());
00239       for (int i=0; i<rows();++i) indices().coeffRef(i) = lhs.indices().coeff(rhs.indices().coeff(i));
00240     }
00241 #endif
00242 
00243   public:
00244 
00249     template<typename Other>
00250     inline PlainPermutationType operator*(const PermutationBase<Other>& other) const
00251     { return PlainPermutationType(internal::PermPermProduct, derived(), other.derived()); }
00252 
00257     template<typename Other>
00258     inline PlainPermutationType operator*(const Transpose<PermutationBase<Other> >& other) const
00259     { return PlainPermutationType(internal::PermPermProduct, *this, other.eval()); }
00260 
00265     template<typename Other> friend
00266     inline PlainPermutationType operator*(const Transpose<PermutationBase<Other> >& other, const PermutationBase& perm)
00267     { return PlainPermutationType(internal::PermPermProduct, other.eval(), perm); }
00268 
00269   protected:
00270 
00271 };
00272 
00287 namespace internal {
00288 template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename IndexType>
00289 struct traits<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime, IndexType> >
00290  : traits<Matrix<IndexType,SizeAtCompileTime,SizeAtCompileTime,0,MaxSizeAtCompileTime,MaxSizeAtCompileTime> >
00291 {
00292   typedef IndexType Index;
00293   typedef Matrix<IndexType, SizeAtCompileTime, 1, 0, MaxSizeAtCompileTime, 1> IndicesType;
00294 };
00295 }
00296 
00297 template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename IndexType>
00298 class PermutationMatrix : public PermutationBase<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime, IndexType> >
00299 {
00300     typedef PermutationBase<PermutationMatrix> Base;
00301     typedef internal::traits<PermutationMatrix> Traits;
00302   public:
00303 
00304     #ifndef EIGEN_PARSED_BY_DOXYGEN
00305     typedef typename Traits::IndicesType IndicesType;
00306     #endif
00307 
00308     inline PermutationMatrix()
00309     {}
00310 
00313     inline PermutationMatrix(int size) : m_indices(size)
00314     {}
00315 
00317     template<typename OtherDerived>
00318     inline PermutationMatrix(const PermutationBase<OtherDerived>& other)
00319       : m_indices(other.indices()) {}
00320 
00321     #ifndef EIGEN_PARSED_BY_DOXYGEN
00322 
00324     inline PermutationMatrix(const PermutationMatrix& other) : m_indices(other.indices()) {}
00325     #endif
00326 
00334     template<typename Other>
00335     explicit inline PermutationMatrix(const MatrixBase<Other>& indices) : m_indices(indices)
00336     {}
00337 
00339     template<typename Other>
00340     explicit PermutationMatrix(const TranspositionsBase<Other>& tr)
00341       : m_indices(tr.size())
00342     {
00343       *this = tr;
00344     }
00345 
00347     template<typename Other>
00348     PermutationMatrix& operator=(const PermutationBase<Other>& other)
00349     {
00350       m_indices = other.indices();
00351       return *this;
00352     }
00353 
00355     template<typename Other>
00356     PermutationMatrix& operator=(const TranspositionsBase<Other>& tr)
00357     {
00358       return Base::operator=(tr.derived());
00359     }
00360 
00361     #ifndef EIGEN_PARSED_BY_DOXYGEN
00362 
00365     PermutationMatrix& operator=(const PermutationMatrix& other)
00366     {
00367       m_indices = other.m_indices;
00368       return *this;
00369     }
00370     #endif
00371 
00373     const IndicesType& indices() const { return m_indices; }
00375     IndicesType& indices() { return m_indices; }
00376 
00377 
00378     /**** multiplication helpers to hopefully get RVO ****/
00379 
00380 #ifndef EIGEN_PARSED_BY_DOXYGEN
00381     template<typename Other>
00382     PermutationMatrix(const Transpose<PermutationBase<Other> >& other)
00383       : m_indices(other.nestedPermutation().size())
00384     {
00385       for (int i=0; i<m_indices.size();++i) m_indices.coeffRef(other.nestedPermutation().indices().coeff(i)) = i;
00386     }
00387     template<typename Lhs,typename Rhs>
00388     PermutationMatrix(internal::PermPermProduct_t, const Lhs& lhs, const Rhs& rhs)
00389       : m_indices(lhs.indices().size())
00390     {
00391       Base::assignProduct(lhs,rhs);
00392     }
00393 #endif
00394 
00395   protected:
00396 
00397     IndicesType m_indices;
00398 };
00399 
00400 
00401 namespace internal {
00402 template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename IndexType, int _PacketAccess>
00403 struct traits<Map<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime, IndexType>,_PacketAccess> >
00404  : traits<Matrix<IndexType,SizeAtCompileTime,SizeAtCompileTime,0,MaxSizeAtCompileTime,MaxSizeAtCompileTime> >
00405 {
00406   typedef IndexType Index;
00407   typedef Map<const Matrix<IndexType, SizeAtCompileTime, 1, 0, MaxSizeAtCompileTime, 1>, _PacketAccess> IndicesType;
00408 };
00409 }
00410 
00411 template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename IndexType, int _PacketAccess>
00412 class Map<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime, IndexType>,_PacketAccess>
00413   : public PermutationBase<Map<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime, IndexType>,_PacketAccess> >
00414 {
00415     typedef PermutationBase<Map> Base;
00416     typedef internal::traits<Map> Traits;
00417   public:
00418 
00419     #ifndef EIGEN_PARSED_BY_DOXYGEN
00420     typedef typename Traits::IndicesType IndicesType;
00421     typedef typename IndicesType::Scalar Index;
00422     #endif
00423 
00424     inline Map(const Index* indices)
00425       : m_indices(indices)
00426     {}
00427 
00428     inline Map(const Index* indices, Index size)
00429       : m_indices(indices,size)
00430     {}
00431 
00433     template<typename Other>
00434     Map& operator=(const PermutationBase<Other>& other)
00435     { return Base::operator=(other.derived()); }
00436 
00438     template<typename Other>
00439     Map& operator=(const TranspositionsBase<Other>& tr)
00440     { return Base::operator=(tr.derived()); }
00441 
00442     #ifndef EIGEN_PARSED_BY_DOXYGEN
00443 
00446     Map& operator=(const Map& other)
00447     {
00448       m_indices = other.m_indices;
00449       return *this;
00450     }
00451     #endif
00452 
00454     const IndicesType& indices() const { return m_indices; }
00456     IndicesType& indices() { return m_indices; }
00457 
00458   protected:
00459 
00460     IndicesType m_indices;
00461 };
00462 
00475 struct PermutationStorage {};
00476 
00477 template<typename _IndicesType> class TranspositionsWrapper;
00478 namespace internal {
00479 template<typename _IndicesType>
00480 struct traits<PermutationWrapper<_IndicesType> >
00481 {
00482   typedef PermutationStorage StorageKind;
00483   typedef typename _IndicesType::Scalar Scalar;
00484   typedef typename _IndicesType::Scalar Index;
00485   typedef _IndicesType IndicesType;
00486   enum {
00487     RowsAtCompileTime = _IndicesType::SizeAtCompileTime,
00488     ColsAtCompileTime = _IndicesType::SizeAtCompileTime,
00489     MaxRowsAtCompileTime = IndicesType::MaxRowsAtCompileTime,
00490     MaxColsAtCompileTime = IndicesType::MaxColsAtCompileTime,
00491     Flags = 0,
00492     CoeffReadCost = _IndicesType::CoeffReadCost
00493   };
00494 };
00495 }
00496 
00497 template<typename _IndicesType>
00498 class PermutationWrapper : public PermutationBase<PermutationWrapper<_IndicesType> >
00499 {
00500     typedef PermutationBase<PermutationWrapper> Base;
00501     typedef internal::traits<PermutationWrapper> Traits;
00502   public:
00503 
00504     #ifndef EIGEN_PARSED_BY_DOXYGEN
00505     typedef typename Traits::IndicesType IndicesType;
00506     #endif
00507 
00508     inline PermutationWrapper(const IndicesType& indices)
00509       : m_indices(indices)
00510     {}
00511 
00513     const typename internal::remove_all<typename IndicesType::Nested>::type&
00514     indices() const { return m_indices; }
00515 
00516   protected:
00517 
00518     typename IndicesType::Nested m_indices;
00519 };
00520 
00523 template<typename Derived, typename PermutationDerived>
00524 inline const internal::permut_matrix_product_retval<PermutationDerived, Derived, OnTheRight>
00525 operator*(const MatrixBase<Derived>& matrix,
00526           const PermutationBase<PermutationDerived> &permutation)
00527 {
00528   return internal::permut_matrix_product_retval
00529            <PermutationDerived, Derived, OnTheRight>
00530            (permutation.derived(), matrix.derived());
00531 }
00532 
00535 template<typename Derived, typename PermutationDerived>
00536 inline const internal::permut_matrix_product_retval
00537                <PermutationDerived, Derived, OnTheLeft>
00538 operator*(const PermutationBase<PermutationDerived> &permutation,
00539           const MatrixBase<Derived>& matrix)
00540 {
00541   return internal::permut_matrix_product_retval
00542            <PermutationDerived, Derived, OnTheLeft>
00543            (permutation.derived(), matrix.derived());
00544 }
00545 
00546 namespace internal {
00547 
00548 template<typename PermutationType, typename MatrixType, int Side, bool Transposed>
00549 struct traits<permut_matrix_product_retval<PermutationType, MatrixType, Side, Transposed> >
00550 {
00551   typedef typename MatrixType::PlainObject ReturnType;
00552 };
00553 
00554 template<typename PermutationType, typename MatrixType, int Side, bool Transposed>
00555 struct permut_matrix_product_retval
00556  : public ReturnByValue<permut_matrix_product_retval<PermutationType, MatrixType, Side, Transposed> >
00557 {
00558     typedef typename remove_all<typename MatrixType::Nested>::type MatrixTypeNestedCleaned;
00559 
00560     permut_matrix_product_retval(const PermutationType& perm, const MatrixType& matrix)
00561       : m_permutation(perm), m_matrix(matrix)
00562     {}
00563 
00564     inline int rows() const { return m_matrix.rows(); }
00565     inline int cols() const { return m_matrix.cols(); }
00566 
00567     template<typename Dest> inline void evalTo(Dest& dst) const
00568     {
00569       const int n = Side==OnTheLeft ? rows() : cols();
00570 
00571       if(is_same<MatrixTypeNestedCleaned,Dest>::value && extract_data(dst) == extract_data(m_matrix))
00572       {
00573         // apply the permutation inplace
00574         Matrix<bool,PermutationType::RowsAtCompileTime,1,0,PermutationType::MaxRowsAtCompileTime> mask(m_permutation.size());
00575         mask.fill(false);
00576         int r = 0;
00577         while(r < m_permutation.size())
00578         {
00579           // search for the next seed
00580           while(r<m_permutation.size() && mask[r]) r++;
00581           if(r>=m_permutation.size())
00582             break;
00583           // we got one, let's follow it until we are back to the seed
00584           int k0 = r++;
00585           int kPrev = k0;
00586           mask.coeffRef(k0) = true;
00587           for(int k=m_permutation.indices().coeff(k0); k!=k0; k=m_permutation.indices().coeff(k))
00588           {
00589                   Block<Dest, Side==OnTheLeft ? 1 : Dest::RowsAtCompileTime, Side==OnTheRight ? 1 : Dest::ColsAtCompileTime>(dst, k)
00590             .swap(Block<Dest, Side==OnTheLeft ? 1 : Dest::RowsAtCompileTime, Side==OnTheRight ? 1 : Dest::ColsAtCompileTime>
00591                        (dst,((Side==OnTheLeft) ^ Transposed) ? k0 : kPrev));
00592 
00593             mask.coeffRef(k) = true;
00594             kPrev = k;
00595           }
00596         }
00597       }
00598       else
00599       {
00600         for(int i = 0; i < n; ++i)
00601         {
00602           Block<Dest, Side==OnTheLeft ? 1 : Dest::RowsAtCompileTime, Side==OnTheRight ? 1 : Dest::ColsAtCompileTime>
00603                (dst, ((Side==OnTheLeft) ^ Transposed) ? m_permutation.indices().coeff(i) : i)
00604 
00605           =
00606 
00607           Block<const MatrixTypeNestedCleaned,Side==OnTheLeft ? 1 : MatrixType::RowsAtCompileTime,Side==OnTheRight ? 1 : MatrixType::ColsAtCompileTime>
00608                (m_matrix, ((Side==OnTheRight) ^ Transposed) ? m_permutation.indices().coeff(i) : i);
00609         }
00610       }
00611     }
00612 
00613   protected:
00614     const PermutationType& m_permutation;
00615     typename MatrixType::Nested m_matrix;
00616 };
00617 
00618 /* Template partial specialization for transposed/inverse permutations */
00619 
00620 template<typename Derived>
00621 struct traits<Transpose<PermutationBase<Derived> > >
00622  : traits<Derived>
00623 {};
00624 
00625 } // end namespace internal
00626 
00627 template<typename Derived>
00628 class Transpose<PermutationBase<Derived> >
00629   : public EigenBase<Transpose<PermutationBase<Derived> > >
00630 {
00631     typedef Derived PermutationType;
00632     typedef typename PermutationType::IndicesType IndicesType;
00633     typedef typename PermutationType::PlainPermutationType PlainPermutationType;
00634   public:
00635 
00636     #ifndef EIGEN_PARSED_BY_DOXYGEN
00637     typedef internal::traits<PermutationType> Traits;
00638     typedef typename Derived::DenseMatrixType DenseMatrixType;
00639     enum {
00640       Flags = Traits::Flags,
00641       CoeffReadCost = Traits::CoeffReadCost,
00642       RowsAtCompileTime = Traits::RowsAtCompileTime,
00643       ColsAtCompileTime = Traits::ColsAtCompileTime,
00644       MaxRowsAtCompileTime = Traits::MaxRowsAtCompileTime,
00645       MaxColsAtCompileTime = Traits::MaxColsAtCompileTime
00646     };
00647     typedef typename Traits::Scalar Scalar;
00648     #endif
00649 
00650     Transpose(const PermutationType& p) : m_permutation(p) {}
00651 
00652     inline int rows() const { return m_permutation.rows(); }
00653     inline int cols() const { return m_permutation.cols(); }
00654 
00655     #ifndef EIGEN_PARSED_BY_DOXYGEN
00656     template<typename DenseDerived>
00657     void evalTo(MatrixBase<DenseDerived>& other) const
00658     {
00659       other.setZero();
00660       for (int i=0; i<rows();++i)
00661         other.coeffRef(i, m_permutation.indices().coeff(i)) = typename DenseDerived::Scalar(1);
00662     }
00663     #endif
00664 
00666     PlainPermutationType eval() const { return *this; }
00667 
00668     DenseMatrixType toDenseMatrix() const { return *this; }
00669 
00672     template<typename OtherDerived> friend
00673     inline const internal::permut_matrix_product_retval<PermutationType, OtherDerived, OnTheRight, true>
00674     operator*(const MatrixBase<OtherDerived>& matrix, const Transpose& trPerm)
00675     {
00676       return internal::permut_matrix_product_retval<PermutationType, OtherDerived, OnTheRight, true>(trPerm.m_permutation, matrix.derived());
00677     }
00678 
00681     template<typename OtherDerived>
00682     inline const internal::permut_matrix_product_retval<PermutationType, OtherDerived, OnTheLeft, true>
00683     operator*(const MatrixBase<OtherDerived>& matrix) const
00684     {
00685       return internal::permut_matrix_product_retval<PermutationType, OtherDerived, OnTheLeft, true>(m_permutation, matrix.derived());
00686     }
00687 
00688     const PermutationType& nestedPermutation() const { return m_permutation; }
00689 
00690   protected:
00691     const PermutationType& m_permutation;
00692 };
00693 
00694 template<typename Derived>
00695 const PermutationWrapper<const Derived> MatrixBase<Derived>::asPermutation() const
00696 {
00697   return derived();
00698 }
00699 
00700 } // end namespace Eigen
00701 
00702 #endif // EIGEN_PERMUTATIONMATRIX_H