[eigen] [PATCH] Toeplitz matrix specialization

[Ricard Marxer Piñón <email@xxxxxxxxxxxxxxxx>]

Here is a first draft of a Toeplitz Matrix specialization.
( http://en.wikipedia.org/wiki/Toeplitz_matrix )

A couple of things:
  - I added a few methods to MatrixBase to easily construct Toeplitz matrices the same way as Diagonal matrices (asToeplitz(), isToeplitz())
  - I added a flag ConstantDiagsBit (maybe it would be better to call it ToeplitzBit, I don't know)

Let me know if I'm going on the right direction and what I did wrong.

Toeplitz matrices allow for specialized solver Levinson-Durbin recursion which exploits its structure.

-- 
ricard
http://www.ricardmarxer.com
http://www.caligraft.com
Index: doc/snippets/MatrixBase_asToeplitz.cpp
===================================================================
--- doc/snippets/MatrixBase_asToeplitz.cpp	(revision 0)
+++ doc/snippets/MatrixBase_asToeplitz.cpp	(revision 0)
@@ -0,0 +1 @@
+cout << Vector3i(2,5,6).asToeplitz() << endl;
Index: doc/snippets/MatrixBase_isToeplitz.cpp
===================================================================
--- doc/snippets/MatrixBase_isToeplitz.cpp	(revision 0)
+++ doc/snippets/MatrixBase_isToeplitz.cpp	(revision 0)
@@ -0,0 +1,6 @@
+Matrix3d m;
+m << 1000, 4000, 5000, 2000, 1000, 4000, 3000, 2002, 1000;
+cout << "Here's the matrix m:" << endl << m << endl;
+cout << "m.isToeplitz() returns: " << m.isToeplitz() << endl;
+cout << "m.isToeplitz(1e-3) returns: " << m.isToeplitz(1e-3) << endl;
+
Index: doc/snippets/MatrixBase_select.cpp
===================================================================
--- doc/snippets/MatrixBase_select.cpp	(revision 0)
+++ doc/snippets/MatrixBase_select.cpp	(revision 0)
@@ -0,0 +1,6 @@
+MatrixXi m(3, 3);
+m << 1, 2, 3, 
+     4, 5, 6, 
+     7, 8, 9;
+m = (m.cwise() >= 5).select(-m, m);
+cout << m << endl;
Index: Eigen/src/Core/ToeplitzMatrix.h
===================================================================
--- Eigen/src/Core/ToeplitzMatrix.h	(revision 0)
+++ Eigen/src/Core/ToeplitzMatrix.h	(revision 0)
@@ -0,0 +1,145 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra. Eigen itself is part of the KDE project.
+//
+// Copyright (C) 2006-2008 Ricard Marxer <email@xxxxxxxxxxxxxxxx>
+//
+// Eigen is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 3 of the License, or (at your option) any later version.
+//
+// Alternatively, you can redistribute it and/or
+// modify it under the terms of the GNU General Public License as
+// published by the Free Software Foundation; either version 2 of
+// the License, or (at your option) any later version.
+//
+// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License and a copy of the GNU General Public License along with
+// Eigen. If not, see <http://www.gnu.org/licenses/>.
+
+#ifndef EIGEN_TOEPLITZMATRIX_H
+#define EIGEN_TOEPLITZMATRIX_H
+
+/** \class ToeplitzMatrix
+  * \nonstableyet 
+  *
+  * \brief Expression of a toeplitz matrix
+  *
+  * \param CoeffsVectorType the type of the vector of Toeplitz coefficients
+  *
+  * This class is an expression of a Toeplitz matrix given the vector of first row's
+  * coefficients and the vector of first col's coefficients. It is the return
+  * type of MatrixBase::toeplitz(const OtherDerived&) and most of the time this is
+  * the only way it is used.
+  *
+  * \sa MatrixBase::toeplitz(const OtherDerived&)
+  */
+template<typename CoeffsVectorType>
+struct ei_traits<ToeplitzMatrix<CoeffsVectorType> >
+{
+  typedef typename CoeffsVectorType::Scalar Scalar;
+  typedef typename ei_nested<CoeffsVectorType>::type CoeffsVectorTypeNested;
+  typedef typename ei_unref<CoeffsVectorTypeNested>::type _CoeffsVectorTypeNested;
+  enum {
+    RowsAtCompileTime = CoeffsVectorType::SizeAtCompileTime,
+    ColsAtCompileTime = CoeffsVectorType::SizeAtCompileTime,
+    MaxRowsAtCompileTime = CoeffsVectorType::MaxSizeAtCompileTime,
+    MaxColsAtCompileTime = CoeffsVectorType::MaxSizeAtCompileTime,
+    Flags = (_CoeffsVectorTypeNested::Flags & HereditaryBits) | Toeplitz,
+    CoeffReadCost = _CoeffsVectorTypeNested::CoeffReadCost
+  };
+};
+
+template<typename CoeffsVectorType>
+class ToeplitzMatrix : ei_no_assignment_operator,
+   public MatrixBase<ToeplitzMatrix<CoeffsVectorType> >
+{
+  public:
+
+    EIGEN_GENERIC_PUBLIC_INTERFACE(ToeplitzMatrix)
+
+    // needed to evaluate a ToeplitzMatrix<Xpr> to a ToeplitzMatrix<NestByValue<Vector> >
+    template<typename OtherCoeffsVectorType>
+    inline ToeplitzMatrix(const ToeplitzMatrix<OtherCoeffsVectorType>& other) : m_first_row(other.toeplitzFirstRow()), m_first_col(other.toeplitzFirstCol())
+    {
+      EIGEN_STATIC_ASSERT_VECTOR_ONLY(CoeffsVectorType);
+      EIGEN_STATIC_ASSERT_VECTOR_ONLY(OtherCoeffsVectorType);
+      ei_assert(m_first_row.size() > 0 && m_first_col.size() > 0);
+    }
+  
+    inline ToeplitzMatrix(const CoeffsVectorType& firstRow, const CoeffsVectorType& firstCol) : m_first_row(firstRow), m_first_col(firstCol)
+    {
+      EIGEN_STATIC_ASSERT_VECTOR_ONLY(CoeffsVectorType);
+      ei_assert(m_first_row.size() > 0 && m_first_col.size() > 0);
+    }
+
+    inline ToeplitzMatrix(const CoeffsVectorType& firstRow) : m_first_row(firstRow), m_first_col(firstRow)
+    {
+      EIGEN_STATIC_ASSERT_VECTOR_ONLY(CoeffsVectorType);
+      ei_assert(m_first_row.size() > 0 && m_first_col.size() > 0);
+    }
+
+
+    inline int rows() const { return m_first_col.size(); }
+    inline int cols() const { return m_first_row.size(); }
+
+    inline const Scalar coeff(int row, int col) const
+    {
+      return row < col ? m_first_row.coeff(col - row) : m_first_col.coeff(row - col);
+    }
+
+    inline const CoeffsVectorType& toeplitzFirstRow() const { return m_first_row; }
+    inline const CoeffsVectorType& toeplitzFirstCol() const { return m_first_col; }
+
+  protected:
+    const typename CoeffsVectorType::Nested m_first_row;
+    const typename CoeffsVectorType::Nested m_first_col;
+};
+
+/** \nonstableyet 
+  * \returns an expression of a toeplitz matrix with *this as vector of toeplitz coefficients
+  *
+  * \only_for_vectors
+  *
+  * \addexample AsToeplitzExample \label How to build a toeplitz matrix from a vector
+  *
+  * Example: \include MatrixBase_asToeplitz.cpp
+  * Output: \verbinclude MatrixBase_asToeplitz.out
+  *
+  * \sa class ToeplitzMatrix, isToeplitz()
+  **/
+template<typename Derived>
+inline const ToeplitzMatrix<Derived>
+MatrixBase<Derived>::asToeplitz() const
+{
+  return derived();
+}
+
+/** \nonstableyet 
+  * \returns true if *this is approximately equal to a toeplitz matrix,
+  *          within the precision given by \a prec.
+  *
+  * Example: \include MatrixBase_isToeplitz.cpp
+  * Output: \verbinclude MatrixBase_isToeplitz.out
+  *
+  * \sa asToeplitz()
+  */
+template<typename Derived>
+bool MatrixBase<Derived>::isToeplitz
+(RealScalar prec) const
+{
+  for(int j = 0; j < cols(); ++j)
+    for(int i = 0; i < j; ++i)
+    {
+      if(!ei_isApprox(coeff(i, j), i < j ? coeff(0, j - i) : coeff(i - j, 0), prec)) return false;
+      if(!ei_isApprox(coeff(j, i), j < i ? coeff(0, i - j) : coeff(j - i, 0), prec)) return false;
+    }
+  return true;
+}
+
+#endif // EIGEN_TOEPLITZMATRIX_H
Index: Eigen/src/Core/MatrixBase.h
===================================================================
--- Eigen/src/Core/MatrixBase.h	(revision 917386)
+++ Eigen/src/Core/MatrixBase.h	(working copy)
@@ -445,6 +445,8 @@
 
     const DiagonalMatrix<Derived> asDiagonal() const;
 
+    const ToeplitzMatrix<Derived> asToeplitz() const;
+
     void fill(const Scalar& value);
     Derived& setConstant(const Scalar& value);
     Derived& setZero();
@@ -476,6 +478,8 @@
                       RealScalar prec = precision<Scalar>()) const;
     bool isUnitary(RealScalar prec = precision<Scalar>()) const;
 
+    bool isToeplitz(RealScalar prec = precision<Scalar>()) const;
+
     template<typename OtherDerived>
     inline bool operator==(const MatrixBase<OtherDerived>& other) const
     { return (cwise() == other).all(); }
Index: Eigen/src/Core/util/Constants.h
===================================================================
--- Eigen/src/Core/util/Constants.h	(revision 917386)
+++ Eigen/src/Core/util/Constants.h	(working copy)
@@ -178,6 +178,11 @@
   * means the expression includes sparse matrices and the sparse path has to be taken. */
 const unsigned int SparseBit = 0x1000;
 
+/** \ingroup flags
+  *
+  * means the all diagonals are constant */
+const unsigned int ConstantDiagsBit = 0x2000;
+
 // list of flags that are inherited by default
 const unsigned int HereditaryBits = RowMajorBit
                                   | EvalBeforeNestingBit
@@ -195,6 +200,7 @@
 const unsigned int UnitUpperTriangular = UpperTriangularBit | UnitDiagBit;
 const unsigned int UnitLowerTriangular = LowerTriangularBit | UnitDiagBit;
 const unsigned int Diagonal = UpperTriangular | LowerTriangular;
+const unsigned int Toeplitz = ConstantDiagsBit;
 
 enum { Aligned, Unaligned };
 enum { ForceAligned, AsRequested };
Index: Eigen/src/Core/util/ForwardDeclarations.h
===================================================================
--- Eigen/src/Core/util/ForwardDeclarations.h	(revision 917386)
+++ Eigen/src/Core/util/ForwardDeclarations.h	(working copy)
@@ -46,6 +46,7 @@
 template<typename BinaryOp,  typename Lhs, typename Rhs>  class CwiseBinaryOp;
 template<typename Lhs, typename Rhs, int ProductMode> class Product;
 template<typename CoeffsVectorType> class DiagonalMatrix;
+template<typename CoeffsVectorType> class ToeplitzMatrix;
 template<typename MatrixType> class DiagonalCoeffs;
 template<typename MatrixType, int PacketAccess = AsRequested> class Map;
 template<typename MatrixType, unsigned int Mode> class Part;
Index: Eigen/Core
===================================================================
--- Eigen/Core	(revision 917386)
+++ Eigen/Core	(working copy)
@@ -147,6 +147,7 @@
 #include "src/Core/CommaInitializer.h"
 #include "src/Core/Part.h"
 #include "src/Core/CacheFriendlyProduct.h"
+#include "src/Core/ToeplitzMatrix.h"
 
 } // namespace Eigen