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

Question

我正在进行大量的科学编程，并使用 Boost.Units 和 Eigen 获得了很好的体验，Boost.Units 提供了数量的编译时维度分析（即用单位标记数量，从而通过经典的物理维度分析捕获了许多错误） 2 线性代数。

然而，Eigen 没有单位的概念，虽然您可以在 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 和模板专业化肯定是可能的。

Answer 1

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

由OP编辑：作为参考，这里是我测试 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;
}

Answer 2

你应该检查这个维基页面：
http://eigen.tuxfamily.org/dox-devel/TopicCustomizingEigen.html

使用原始数据类型之外的 Eigen 需要做一些工作，但这通常是可能的。

Answer 3

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

例如，对于与 * 乘法运算符一起使用的 Boost 单位自定义类型，对于任意 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 * 运算符将返回与输入相同的“类型” - 这就是问题所在。

另一种选择是将特征矩阵嵌入到数量中，而不是将数量嵌入到矩阵中。

Answer 4

具有统一单位的向量/矩阵的简单情况可以通过调整特征矩阵的值类型来部分解决，但这需要大量的用户以专门特征的形式进行干预，甚至这还不足以满足所有操作。

但是，如果您想解决包含非均匀单位的矩阵的一般情况（即每个条目可以具有 affigeren 单位），那么通过自定义矩阵的值类型将不起作用。

您需要一种像此处介绍的方法，该方法位于线性代数库之上并提供单位注释： https://m.youtube.com/watch?v=4LmMwhM8ODI