Eigen  3.3.3
PermutationMatrix.h
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-2015 Gael Guennebaud <gael.guennebaud@inria.fr>
00006 //
00007 // This Source Code Form is subject to the terms of the Mozilla
00008 // Public License v. 2.0. If a copy of the MPL was not distributed
00009 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
00010 
00011 #ifndef EIGEN_PERMUTATIONMATRIX_H
00012 #define EIGEN_PERMUTATIONMATRIX_H
00013 
00014 namespace Eigen { 
00015 
00016 namespace internal {
00017 
00018 enum PermPermProduct_t {PermPermProduct};
00019 
00020 } // end namespace internal
00021 
00045 template<typename Derived>
00046 class PermutationBase : public EigenBase<Derived>
00047 {
00048     typedef internal::traits<Derived> Traits;
00049     typedef EigenBase<Derived> Base;
00050   public:
00051 
00052     #ifndef EIGEN_PARSED_BY_DOXYGEN
00053     typedef typename Traits::IndicesType IndicesType;
00054     enum {
00055       Flags = Traits::Flags,
00056       RowsAtCompileTime = Traits::RowsAtCompileTime,
00057       ColsAtCompileTime = Traits::ColsAtCompileTime,
00058       MaxRowsAtCompileTime = Traits::MaxRowsAtCompileTime,
00059       MaxColsAtCompileTime = Traits::MaxColsAtCompileTime
00060     };
00061     typedef typename Traits::StorageIndex StorageIndex;
00062     typedef Matrix<StorageIndex,RowsAtCompileTime,ColsAtCompileTime,0,MaxRowsAtCompileTime,MaxColsAtCompileTime>
00063             DenseMatrixType;
00064     typedef PermutationMatrix<IndicesType::SizeAtCompileTime,IndicesType::MaxSizeAtCompileTime,StorageIndex>
00065             PlainPermutationType;
00066     typedef PlainPermutationType PlainObject;
00067     using Base::derived;
00068     typedef Inverse<Derived> InverseReturnType;
00069     typedef void Scalar;
00070     #endif
00071 
00073     template<typename OtherDerived>
00074     Derived& operator=(const PermutationBase<OtherDerived>& other)
00075     {
00076       indices() = other.indices();
00077       return derived();
00078     }
00079 
00081     template<typename OtherDerived>
00082     Derived& operator=(const TranspositionsBase<OtherDerived>& tr)
00083     {
00084       setIdentity(tr.size());
00085       for(Index k=size()-1; k>=0; --k)
00086         applyTranspositionOnTheRight(k,tr.coeff(k));
00087       return derived();
00088     }
00089 
00090     #ifndef EIGEN_PARSED_BY_DOXYGEN
00091 
00094     Derived& operator=(const PermutationBase& other)
00095     {
00096       indices() = other.indices();
00097       return derived();
00098     }
00099     #endif
00100 
00102     inline Index rows() const { return Index(indices().size()); }
00103 
00105     inline Index cols() const { return Index(indices().size()); }
00106 
00108     inline Index size() const { return Index(indices().size()); }
00109 
00110     #ifndef EIGEN_PARSED_BY_DOXYGEN
00111     template<typename DenseDerived>
00112     void evalTo(MatrixBase<DenseDerived>& other) const
00113     {
00114       other.setZero();
00115       for (Index i=0; i<rows(); ++i)
00116         other.coeffRef(indices().coeff(i),i) = typename DenseDerived::Scalar(1);
00117     }
00118     #endif
00119 
00124     DenseMatrixType toDenseMatrix() const
00125     {
00126       return derived();
00127     }
00128 
00130     const IndicesType& indices() const { return derived().indices(); }
00132     IndicesType& indices() { return derived().indices(); }
00133 
00136     inline void resize(Index newSize)
00137     {
00138       indices().resize(newSize);
00139     }
00140 
00142     void setIdentity()
00143     {
00144       StorageIndex n = StorageIndex(size());
00145       for(StorageIndex i = 0; i < n; ++i)
00146         indices().coeffRef(i) = i;
00147     }
00148 
00151     void setIdentity(Index newSize)
00152     {
00153       resize(newSize);
00154       setIdentity();
00155     }
00156 
00166     Derived& applyTranspositionOnTheLeft(Index i, Index j)
00167     {
00168       eigen_assert(i>=0 && j>=0 && i<size() && j<size());
00169       for(Index k = 0; k < size(); ++k)
00170       {
00171         if(indices().coeff(k) == i) indices().coeffRef(k) = StorageIndex(j);
00172         else if(indices().coeff(k) == j) indices().coeffRef(k) = StorageIndex(i);
00173       }
00174       return derived();
00175     }
00176 
00185     Derived& applyTranspositionOnTheRight(Index i, Index j)
00186     {
00187       eigen_assert(i>=0 && j>=0 && i<size() && j<size());
00188       std::swap(indices().coeffRef(i), indices().coeffRef(j));
00189       return derived();
00190     }
00191 
00196     inline InverseReturnType inverse() const
00197     { return InverseReturnType(derived()); }
00202     inline InverseReturnType transpose() const
00203     { return InverseReturnType(derived()); }
00204 
00205     /**** multiplication helpers to hopefully get RVO ****/
00206 
00207   
00208 #ifndef EIGEN_PARSED_BY_DOXYGEN
00209   protected:
00210     template<typename OtherDerived>
00211     void assignTranspose(const PermutationBase<OtherDerived>& other)
00212     {
00213       for (Index i=0; i<rows();++i) indices().coeffRef(other.indices().coeff(i)) = i;
00214     }
00215     template<typename Lhs,typename Rhs>
00216     void assignProduct(const Lhs& lhs, const Rhs& rhs)
00217     {
00218       eigen_assert(lhs.cols() == rhs.rows());
00219       for (Index i=0; i<rows();++i) indices().coeffRef(i) = lhs.indices().coeff(rhs.indices().coeff(i));
00220     }
00221 #endif
00222 
00223   public:
00224 
00229     template<typename Other>
00230     inline PlainPermutationType operator*(const PermutationBase<Other>& other) const
00231     { return PlainPermutationType(internal::PermPermProduct, derived(), other.derived()); }
00232 
00237     template<typename Other>
00238     inline PlainPermutationType operator*(const InverseImpl<Other,PermutationStorage>& other) const
00239     { return PlainPermutationType(internal::PermPermProduct, *this, other.eval()); }
00240 
00245     template<typename Other> friend
00246     inline PlainPermutationType operator*(const InverseImpl<Other, PermutationStorage>& other, const PermutationBase& perm)
00247     { return PlainPermutationType(internal::PermPermProduct, other.eval(), perm); }
00248     
00253     Index determinant() const
00254     {
00255       Index res = 1;
00256       Index n = size();
00257       Matrix<bool,RowsAtCompileTime,1,0,MaxRowsAtCompileTime> mask(n);
00258       mask.fill(false);
00259       Index r = 0;
00260       while(r < n)
00261       {
00262         // search for the next seed
00263         while(r<n && mask[r]) r++;
00264         if(r>=n)
00265           break;
00266         // we got one, let's follow it until we are back to the seed
00267         Index k0 = r++;
00268         mask.coeffRef(k0) = true;
00269         for(Index k=indices().coeff(k0); k!=k0; k=indices().coeff(k))
00270         {
00271           mask.coeffRef(k) = true;
00272           res = -res;
00273         }
00274       }
00275       return res;
00276     }
00277 
00278   protected:
00279 
00280 };
00281 
00282 namespace internal {
00283 template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename _StorageIndex>
00284 struct traits<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime, _StorageIndex> >
00285  : traits<Matrix<_StorageIndex,SizeAtCompileTime,SizeAtCompileTime,0,MaxSizeAtCompileTime,MaxSizeAtCompileTime> >
00286 {
00287   typedef PermutationStorage StorageKind;
00288   typedef Matrix<_StorageIndex, SizeAtCompileTime, 1, 0, MaxSizeAtCompileTime, 1> IndicesType;
00289   typedef _StorageIndex StorageIndex;
00290   typedef void Scalar;
00291 };
00292 }
00293 
00307 template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename _StorageIndex>
00308 class PermutationMatrix : public PermutationBase<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime, _StorageIndex> >
00309 {
00310     typedef PermutationBase<PermutationMatrix> Base;
00311     typedef internal::traits<PermutationMatrix> Traits;
00312   public:
00313 
00314     typedef const PermutationMatrix& Nested;
00315 
00316     #ifndef EIGEN_PARSED_BY_DOXYGEN
00317     typedef typename Traits::IndicesType IndicesType;
00318     typedef typename Traits::StorageIndex StorageIndex;
00319     #endif
00320 
00321     inline PermutationMatrix()
00322     {}
00323 
00326     explicit inline PermutationMatrix(Index size) : m_indices(size)
00327     {
00328       eigen_internal_assert(size <= NumTraits<StorageIndex>::highest());
00329     }
00330 
00332     template<typename OtherDerived>
00333     inline PermutationMatrix(const PermutationBase<OtherDerived>& other)
00334       : m_indices(other.indices()) {}
00335 
00336     #ifndef EIGEN_PARSED_BY_DOXYGEN
00337 
00339     inline PermutationMatrix(const PermutationMatrix& other) : m_indices(other.indices()) {}
00340     #endif
00341 
00349     template<typename Other>
00350     explicit inline PermutationMatrix(const MatrixBase<Other>& indices) : m_indices(indices)
00351     {}
00352 
00354     template<typename Other>
00355     explicit PermutationMatrix(const TranspositionsBase<Other>& tr)
00356       : m_indices(tr.size())
00357     {
00358       *this = tr;
00359     }
00360 
00362     template<typename Other>
00363     PermutationMatrix& operator=(const PermutationBase<Other>& other)
00364     {
00365       m_indices = other.indices();
00366       return *this;
00367     }
00368 
00370     template<typename Other>
00371     PermutationMatrix& operator=(const TranspositionsBase<Other>& tr)
00372     {
00373       return Base::operator=(tr.derived());
00374     }
00375 
00376     #ifndef EIGEN_PARSED_BY_DOXYGEN
00377 
00380     PermutationMatrix& operator=(const PermutationMatrix& other)
00381     {
00382       m_indices = other.m_indices;
00383       return *this;
00384     }
00385     #endif
00386 
00388     const IndicesType& indices() const { return m_indices; }
00390     IndicesType& indices() { return m_indices; }
00391 
00392 
00393     /**** multiplication helpers to hopefully get RVO ****/
00394 
00395 #ifndef EIGEN_PARSED_BY_DOXYGEN
00396     template<typename Other>
00397     PermutationMatrix(const InverseImpl<Other,PermutationStorage>& other)
00398       : m_indices(other.derived().nestedExpression().size())
00399     {
00400       eigen_internal_assert(m_indices.size() <= NumTraits<StorageIndex>::highest());
00401       StorageIndex end = StorageIndex(m_indices.size());
00402       for (StorageIndex i=0; i<end;++i)
00403         m_indices.coeffRef(other.derived().nestedExpression().indices().coeff(i)) = i;
00404     }
00405     template<typename Lhs,typename Rhs>
00406     PermutationMatrix(internal::PermPermProduct_t, const Lhs& lhs, const Rhs& rhs)
00407       : m_indices(lhs.indices().size())
00408     {
00409       Base::assignProduct(lhs,rhs);
00410     }
00411 #endif
00412 
00413   protected:
00414 
00415     IndicesType m_indices;
00416 };
00417 
00418 
00419 namespace internal {
00420 template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename _StorageIndex, int _PacketAccess>
00421 struct traits<Map<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime, _StorageIndex>,_PacketAccess> >
00422  : traits<Matrix<_StorageIndex,SizeAtCompileTime,SizeAtCompileTime,0,MaxSizeAtCompileTime,MaxSizeAtCompileTime> >
00423 {
00424   typedef PermutationStorage StorageKind;
00425   typedef Map<const Matrix<_StorageIndex, SizeAtCompileTime, 1, 0, MaxSizeAtCompileTime, 1>, _PacketAccess> IndicesType;
00426   typedef _StorageIndex StorageIndex;
00427   typedef void Scalar;
00428 };
00429 }
00430 
00431 template<int SizeAtCompileTime, int MaxSizeAtCompileTime, typename _StorageIndex, int _PacketAccess>
00432 class Map<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime, _StorageIndex>,_PacketAccess>
00433   : public PermutationBase<Map<PermutationMatrix<SizeAtCompileTime, MaxSizeAtCompileTime, _StorageIndex>,_PacketAccess> >
00434 {
00435     typedef PermutationBase<Map> Base;
00436     typedef internal::traits<Map> Traits;
00437   public:
00438 
00439     #ifndef EIGEN_PARSED_BY_DOXYGEN
00440     typedef typename Traits::IndicesType IndicesType;
00441     typedef typename IndicesType::Scalar StorageIndex;
00442     #endif
00443 
00444     inline Map(const StorageIndex* indicesPtr)
00445       : m_indices(indicesPtr)
00446     {}
00447 
00448     inline Map(const StorageIndex* indicesPtr, Index size)
00449       : m_indices(indicesPtr,size)
00450     {}
00451 
00453     template<typename Other>
00454     Map& operator=(const PermutationBase<Other>& other)
00455     { return Base::operator=(other.derived()); }
00456 
00458     template<typename Other>
00459     Map& operator=(const TranspositionsBase<Other>& tr)
00460     { return Base::operator=(tr.derived()); }
00461 
00462     #ifndef EIGEN_PARSED_BY_DOXYGEN
00463 
00466     Map& operator=(const Map& other)
00467     {
00468       m_indices = other.m_indices;
00469       return *this;
00470     }
00471     #endif
00472 
00474     const IndicesType& indices() const { return m_indices; }
00476     IndicesType& indices() { return m_indices; }
00477 
00478   protected:
00479 
00480     IndicesType m_indices;
00481 };
00482 
00483 template<typename _IndicesType> class TranspositionsWrapper;
00484 namespace internal {
00485 template<typename _IndicesType>
00486 struct traits<PermutationWrapper<_IndicesType> >
00487 {
00488   typedef PermutationStorage StorageKind;
00489   typedef void Scalar;
00490   typedef typename _IndicesType::Scalar StorageIndex;
00491   typedef _IndicesType IndicesType;
00492   enum {
00493     RowsAtCompileTime = _IndicesType::SizeAtCompileTime,
00494     ColsAtCompileTime = _IndicesType::SizeAtCompileTime,
00495     MaxRowsAtCompileTime = IndicesType::MaxSizeAtCompileTime,
00496     MaxColsAtCompileTime = IndicesType::MaxSizeAtCompileTime,
00497     Flags = 0
00498   };
00499 };
00500 }
00501 
00513 template<typename _IndicesType>
00514 class PermutationWrapper : public PermutationBase<PermutationWrapper<_IndicesType> >
00515 {
00516     typedef PermutationBase<PermutationWrapper> Base;
00517     typedef internal::traits<PermutationWrapper> Traits;
00518   public:
00519 
00520     #ifndef EIGEN_PARSED_BY_DOXYGEN
00521     typedef typename Traits::IndicesType IndicesType;
00522     #endif
00523 
00524     inline PermutationWrapper(const IndicesType& indices)
00525       : m_indices(indices)
00526     {}
00527 
00529     const typename internal::remove_all<typename IndicesType::Nested>::type&
00530     indices() const { return m_indices; }
00531 
00532   protected:
00533 
00534     typename IndicesType::Nested m_indices;
00535 };
00536 
00537 
00540 template<typename MatrixDerived, typename PermutationDerived>
00541 EIGEN_DEVICE_FUNC
00542 const Product<MatrixDerived, PermutationDerived, AliasFreeProduct>
00543 operator*(const MatrixBase<MatrixDerived> &matrix,
00544           const PermutationBase<PermutationDerived>& permutation)
00545 {
00546   return Product<MatrixDerived, PermutationDerived, AliasFreeProduct>
00547             (matrix.derived(), permutation.derived());
00548 }
00549 
00552 template<typename PermutationDerived, typename MatrixDerived>
00553 EIGEN_DEVICE_FUNC
00554 const Product<PermutationDerived, MatrixDerived, AliasFreeProduct>
00555 operator*(const PermutationBase<PermutationDerived> &permutation,
00556           const MatrixBase<MatrixDerived>& matrix)
00557 {
00558   return Product<PermutationDerived, MatrixDerived, AliasFreeProduct>
00559             (permutation.derived(), matrix.derived());
00560 }
00561 
00562 
00563 template<typename PermutationType>
00564 class InverseImpl<PermutationType, PermutationStorage>
00565   : public EigenBase<Inverse<PermutationType> >
00566 {
00567     typedef typename PermutationType::PlainPermutationType PlainPermutationType;
00568     typedef internal::traits<PermutationType> PermTraits;
00569   protected:
00570     InverseImpl() {}
00571   public:
00572     typedef Inverse<PermutationType> InverseType;
00573     using EigenBase<Inverse<PermutationType> >::derived;
00574 
00575     #ifndef EIGEN_PARSED_BY_DOXYGEN
00576     typedef typename PermutationType::DenseMatrixType DenseMatrixType;
00577     enum {
00578       RowsAtCompileTime = PermTraits::RowsAtCompileTime,
00579       ColsAtCompileTime = PermTraits::ColsAtCompileTime,
00580       MaxRowsAtCompileTime = PermTraits::MaxRowsAtCompileTime,
00581       MaxColsAtCompileTime = PermTraits::MaxColsAtCompileTime
00582     };
00583     #endif
00584 
00585     #ifndef EIGEN_PARSED_BY_DOXYGEN
00586     template<typename DenseDerived>
00587     void evalTo(MatrixBase<DenseDerived>& other) const
00588     {
00589       other.setZero();
00590       for (Index i=0; i<derived().rows();++i)
00591         other.coeffRef(i, derived().nestedExpression().indices().coeff(i)) = typename DenseDerived::Scalar(1);
00592     }
00593     #endif
00594 
00596     PlainPermutationType eval() const { return derived(); }
00597 
00598     DenseMatrixType toDenseMatrix() const { return derived(); }
00599 
00602     template<typename OtherDerived> friend
00603     const Product<OtherDerived, InverseType, AliasFreeProduct>
00604     operator*(const MatrixBase<OtherDerived>& matrix, const InverseType& trPerm)
00605     {
00606       return Product<OtherDerived, InverseType, AliasFreeProduct>(matrix.derived(), trPerm.derived());
00607     }
00608 
00611     template<typename OtherDerived>
00612     const Product<InverseType, OtherDerived, AliasFreeProduct>
00613     operator*(const MatrixBase<OtherDerived>& matrix) const
00614     {
00615       return Product<InverseType, OtherDerived, AliasFreeProduct>(derived(), matrix.derived());
00616     }
00617 };
00618 
00619 template<typename Derived>
00620 const PermutationWrapper<const Derived> MatrixBase<Derived>::asPermutation() const
00621 {
00622   return derived();
00623 }
00624 
00625 namespace internal {
00626 
00627 template<> struct AssignmentKind<DenseShape,PermutationShape> { typedef EigenBase2EigenBase Kind; };
00628 
00629 } // end namespace internal
00630 
00631 } // end namespace Eigen
00632 
00633 #endif // EIGEN_PERMUTATIONMATRIX_H
 All Classes Functions Variables Typedefs Enumerations Enumerator Friends