将线性代数库与Boost::Units结合使用

Question

我在进行一些科学编程，并使用了Boost.Units和Eigen 2，都取得了很好的效果。

Boost.Units可以为量提供编译时的量纲分析（即标记带单位的量并通过经典物理量纲分析来捕获许多错误），而我在线性代数方面则使用Eigen 2。

然而，Eigen没有单位的概念。虽然你可以在矩阵中设置标量量，并将两个量相乘，但它期望相乘后的类型保持不变，而这显然对于单位是不正确的。例如，以下代码：

using boost::units::quantity;
namespace si = boost::units::si;
Eigen::Matrix< quantity< si::length >, 2, 1 > meter_vector;
quantity< si::area > norm = meter_vector.squaredNorm();

尽管它在逻辑上是正确的，但不起作用。

是否有任何支持单位的矩阵库？我知道这在过去实现起来非常困难，C++11和decltype会使这更容易，但使用C++03和模板特化也是可能的。

- thiton

Answer 1

我相信Blitz++支持大部分Boost.Units功能。

编辑者注：以下是我用来测试Blitz矩阵乘法功能的完整代码：

#include <blitz/array.h>
#include <boost/units/systems/si/area.hpp>
#include <boost/units/systems/si/length.hpp>
#include <boost/units/quantity.hpp>

using boost::units::quantity;
namespace si = boost::units::si;

namespace blitz {
template< typename U1, typename T1, typename U2, typename T2>
struct Multiply< quantity<U1,T1>, quantity<U2,T2> >
{
    typedef typename boost::units::multiply_typeof_helper< quantity<U1,T1>, quantity<U2,T2> >::type T_numtype;

    static inline T_numtype apply( quantity<U1,T1> a, quantity<U2,T2> b ) { return a*b; }
};

}

using namespace blitz;

int main() {
    Array< quantity<si::length>, 1 > matrix;
    Array< quantity<si::area>, 1 > area;
    area = matrix * matrix;
    return 0;
}

- Matthias Schabel

---

顺便说一下，因为我自己不得不搜索一下：The blitz manual 3.7.1告诉你如何提升用户定义的类型。感谢提示。 - thiton

Answer 2

你应该查看这个Wiki页面： http://eigen.tuxfamily.org/dox-devel/TopicCustomizingEigen.html Eigen需要一些工作才能使用除原始数据类型之外的其他类型，但通常是可能的。

- mrolf

---

谢谢提示。已经阅读了页面并按照提示操作。关键是，operator+ 运行良好，但例如 operator* 是错误的，因为米 * 米 不是一个米。 - thiton

Answer 3

使用标准Eigen库插件选项的困难在于，现有的操作符+、-、*等需要被替换为Boost Units量才能使用。

例如，对于一个Boost units自定义类型与乘法运算符*一起工作，对于任意CUSTOM_TYPE，它需要像这样：

template<class X,class Y>
CUSTOM_TYPE<typename boost::units::multiply_typeof_helper<X,Y>::type>
operator*(const CUSTOM_TYPE<X>& x,const CUSTOM_TYPE<Y>& y)
{
    typedef typename boost::units::multiply_typeof_helper<X,Y>::type    type;

    return CUSTOM_TYPE<type>( ... );
}

注意返回类型与输入类型不同。在这里，您使用模板助手multiply_typeof_helper创建返回类型。这是因为将米乘以秒不会给您任何单位的数量。但是，默认的Eigen *运算符将返回与输入相同的“类型”-这是问题所在。

另一个选择是将Eigen矩阵嵌入到量中，而不是将量嵌入到矩阵中。

- Deniz

Answer 4

对于具有统一单位的向量/矩阵的简单情况，可以通过调整Eigen矩阵的值类型来部分解决，但这需要用户在特化特性方面进行相当大的干预，而且即使如此，对于所有操作来说仍然不足够。

然而，如果你想解决包含非统一单位（即每个条目都可以有不同单位）的矩阵的一般情况，通过自定义矩阵的值类型是行不通的。

你需要像这里介绍的方法一样，构建在线性代数库之上，并提供单位注释：https://m.youtube.com/watch?v=4LmMwhM8ODI

- Daniel