[Eigen中文文档] 扩展/自定义Eigen（一）

文档总目录

1. 扩展 MatrixBase（包括其他类）

英文原文(Extending MatrixBase (and other classes))

在本节中，将介绍如何向MatrixBase添加自定义方法。由于所有表达式和矩阵类型都继承自MatrixBase，因此将方法添加到MatrixBase会立即使其对所有表达式可用！一个典型的用例是，使Eigen与另一个API兼容。

你肯定知道在C++中无法向现有类添加方法。Eigen里的技巧是在MatrixBase的声明中包含由预处理器宏EIGEN_MATRIXBASE_PLUGIN 定义的文件，如下：

class MatrixBase {
  // ...
  #ifdef EIGEN_MATRIXBASE_PLUGIN
  #include EIGEN_MATRIXBASE_PLUGIN
  #endif
};

因此，要使用自己的方法扩展 MatrixBase，只需创建一个包含方法声明的文件，并在包含任何 Eigen 的头文件之前定义 EIGEN_MATRIXBASE_PLUGIN。

你可以通过定义类似命名的预处理器符号来扩展Eigen中使用的许多其他类。例如，如果要扩展ArrayBase类，则定义EIGEN_ARRAYBASE_PLUGIN。可以在 预处理器指令 中找到可以以此方式扩展的所有类以及相应的预处理器符号的完整列表。

以下是向 MatrixBase 添加方法的扩展文件示例：

MatrixBaseAddons.h

inline Scalar at(uint i, uint j) const { return this->operator()(i,j); }
inline Scalar& at(uint i, uint j) { return this->operator()(i,j); }
inline Scalar at(uint i) const { return this->operator[](i); }
inline Scalar& at(uint i) { return this->operator[](i); }
 
inline RealScalar squaredLength() const { return squaredNorm(); }
inline RealScalar length() const { return norm(); }
inline RealScalar invLength(void) const { return fast_inv_sqrt(squaredNorm()); }
 
template<typename OtherDerived>
inline Scalar squaredDistanceTo(const MatrixBase<OtherDerived>& other) const
{ return (derived() - other.derived()).squaredNorm(); }
 
template<typename OtherDerived>
inline RealScalar distanceTo(const MatrixBase<OtherDerived>& other) const
{ return internal::sqrt(derived().squaredDistanceTo(other)); }
 
inline void scaleTo(RealScalar l) { RealScalar vl = norm(); if (vl>1e-9) derived() *= (l/vl); }
 
inline Transpose<Derived> transposed() {return this->transpose();}
inline const Transpose<Derived> transposed() const {return this->transpose();}
 
inline uint minComponentId(void) const  { int i; this->minCoeff(&i); return i; }
inline uint maxComponentId(void) const  { int i; this->maxCoeff(&i); return i; }
 
template<typename OtherDerived>
void makeFloor(const MatrixBase<OtherDerived>& other) { derived() = derived().cwiseMin(other.derived()); }
template<typename OtherDerived>
void makeCeil(const MatrixBase<OtherDerived>& other) { derived() = derived().cwiseMax(other.derived()); }
 
const CwiseBinaryOp<internal::scalar_sum_op<Scalar>, const Derived, const ConstantReturnType>
operator+(const Scalar& scalar) const
{ return CwiseBinaryOp<internal::scalar_sum_op<Scalar>, const Derived, const ConstantReturnType>(derived(), Constant(rows(),cols(),scalar)); }
 
friend const CwiseBinaryOp<internal::scalar_sum_op<Scalar>, const ConstantReturnType, Derived>
operator+(const Scalar& scalar, const MatrixBase<Derived>& mat)
{ return CwiseBinaryOp<internal::scalar_sum_op<Scalar>, const ConstantReturnType, Derived>(Constant(rows(),cols(),scalar), mat.derived()); }

然后可以在 config.h 或项目的任何先决条件头文件中进行以下声明：

#define EIGEN_MATRIXBASE_PLUGIN "MatrixBaseAddons.h"

2. 继承 Matrix

英文原文(Inheriting from Matrix)

在从Matrix继承之前，请确定使用 EIGEN_MATRIX_PLUGIN 不是您真正想要的（请参见前一节）。如果只需要向Matrix添加几个成员，那么继承Matrix就是正确的方法。

当有多层继承关系，例如 MyVerySpecificVector1, MyVerySpecificVector2 -> MyVector1 -> Matrix 和 MyVerySpecificVector3, MyVerySpecificVector4 -> MyVector2 -> Matrix ，这时需要继承Matrix。

为了使对象在 Eigen 框架内工作，需要在继承的类中定义一些成员，如下：

#include <Eigen/Core>
#include <iostream>
 
class MyVectorType : public Eigen::VectorXd
{
public:
    MyVectorType(void):Eigen::VectorXd() {}
 
    // This constructor allows you to construct MyVectorType from Eigen expressions
    template<typename OtherDerived>
    MyVectorType(const Eigen::MatrixBase<OtherDerived>& other)
        : Eigen::VectorXd(other)
    { }
 
    // This method allows you to assign Eigen expressions to MyVectorType
    template<typename OtherDerived>
    MyVectorType& operator=(const Eigen::MatrixBase <OtherDerived>& other)
    {
        this->Eigen::VectorXd::operator=(other);
        return *this;
    }
};
 
int main()
{
    MyVectorType v = MyVectorType::Ones(4);
    v(2) += 10;
    v = 2 * v;
    std::cout << v.transpose() << std::endl;
}

输出：

2  2 22  2

如果不提供这些方法，可能会遇到这种错误：

error: no match for ‘operator=’ in ‘v = Eigen::operator*(
const Eigen::MatrixBase<Eigen::Matrix<double, -0x000000001, 1, 0, -0x000000001, 1> >::Scalar&, 
const Eigen::MatrixBase<Eigen::Matrix<double, -0x000000001, 1> >::StorageBaseType&)
(((const Eigen::MatrixBase<Eigen::Matrix<double, -0x000000001, 1> >::StorageBaseType&)
((const Eigen::MatrixBase<Eigen::Matrix<double, -0x000000001, 1> >::StorageBaseType*)(& v))))’

3. 使用自定义标量类型

英文原文(Using custom scalar types)

默认情况下，Eigen目前支持标准浮点类型（float、double、std::complex<float>、std::complex<double>、long double），以及所有原生的整型（int、unsigned int、short等）和bool类型。在x86-64系统上，long double允许本地强制使用带扩展精度的x87寄存器（相对于SSE）。

为了添加对自定义类型 T 的支持，需要：

• 确保类型 T 支持常用运算符（+、-、*、/等）
• 添加 struct Eigen::NumTraits<T> 的特殊化（请参阅 NumTraits）
• 定义对自定义类型有意义的数学函数。包括标准的 sqrt、pow、sin、tan、conj、real、imag 等，以及 Eigen 特定的 abs2。(请参阅文件 Eigen/src/Core/MathFunctions.h)

数学函数应该定义在与 T 相同的命名空间中，或者定义在 std 命名空间中，但不建议使用第二种方法。

下面是一个添加对 Adolc 的 adouble 类型支持的具体示例。 Adolc 是一个自动微分库。 adouble 类型基本上是跟踪任意数量的偏导数值的实值。

#ifndef ADOLCSUPPORT_H
#define ADOLCSUPPORT_H
 
#define ADOLC_TAPELESS
#include <adolc/adouble.h>
#include <Eigen/Core>
 
namespace Eigen {
 
template<> struct NumTraits<adtl::adouble>
 : NumTraits<double> // permits to get the epsilon, dummy_precision, lowest, highest functions
{
  typedef adtl::adouble Real;
  typedef adtl::adouble NonInteger;
  typedef adtl::adouble Nested;
 
  enum {
    IsComplex = 0,
    IsInteger = 0,
    IsSigned = 1,
    RequireInitialization = 1,
    ReadCost = 1,
    AddCost = 3,
    MulCost = 3
  };
};
 
}
 
namespace adtl {
 
inline const adouble& conj(const adouble& x)  { return x; }
inline const adouble& real(const adouble& x)  { return x; }
inline adouble imag(const adouble&)    { return 0.; }
inline adouble abs(const adouble&  x)  { return fabs(x); }
inline adouble abs2(const adouble& x)  { return x*x; }
 
}
 
#endif // ADOLCSUPPORT_H

这个例子增加了对GMP库中mpq_class类型的支持。它特别展示了如何在LU分解过程中改变Eigen选择最佳主元的方法。它选择具有最高分数的系数作为主元，在默认情况下分数是数字的绝对值，但我们可以定义不同的分数，比如更喜欢具有更紧凑表示的主元（这只是一个例子，不是建议）。请注意，分数应始终为非负数，并且只允许零具有零的分数。此外，这可能会与不精确标量类型的阈值产生不良交互。

#include <gmpxx.h>
#include <Eigen/Core>
#include <boost/operators.hpp>
 
namespace Eigen {
  template<> struct NumTraits<mpq_class> : GenericNumTraits<mpq_class>
  {
    typedef mpq_class Real;
    typedef mpq_class NonInteger;
    typedef mpq_class Nested;
 
    static inline Real epsilon() { return 0; }
    static inline Real dummy_precision() { return 0; }
    static inline int digits10() { return 0; }
 
    enum {
      IsInteger = 0,
      IsSigned = 1,
      IsComplex = 0,
      RequireInitialization = 1,
      ReadCost = 6,
      AddCost = 150,
      MulCost = 100
    };
  };
 
  namespace internal {
 
    template<> struct scalar_score_coeff_op<mpq_class> {
      struct result_type : boost::totally_ordered1<result_type> {
        std::size_t len;
        result_type(int i = 0) : len(i) {} // Eigen uses Score(0) and Score()
        result_type(mpq_class const& q) :
          len(mpz_size(q.get_num_mpz_t())+
              mpz_size(q.get_den_mpz_t())-1) {}
        friend bool operator<(result_type x, result_type y) {
          // 0 is the worst possible pivot
          if (x.len == 0) return y.len > 0;
          if (y.len == 0) return false;
          // Prefer a pivot with a small representation
          return x.len > y.len;
        }
        friend bool operator==(result_type x, result_type y) {
          // Only used to test if the score is 0
          return x.len == y.len;
        }
      };
      result_type operator()(mpq_class const& x) const { return x; }
    };
  }
}