Combining a linear algebra library with Boost::Units

Question

I'm doing a good amount of scientific programming and made very good experiences with both Boost.Units, which provides compile-time dimensional analysis for quantities (i.e. tags quantities with units and thereby catches many errors with classical physical dimension analysis) and using Eigen 2 for linear algebra.

However, Eigen has no concept of units, and while you can set the scalar quantities in matrices for Eigen, it expects that the multiplication of two quantities yields the same type, which is obviously untrue for units. For example, code like:

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();

does not work, even though it is logically correct.

Is there any matrix library that supports units? I know that this would have been notoriously hard to implement in the past, and C++11 and decltype will make that much easier, but it surely was possible with C++03 and template specializations.

Cannery

Answer 1

I believe Blitz++ supports much of Boost.Units functionality.

Edit by the OP: For the reference here is the full test code with which I tested the Blitz matrix multiplication functionality:

#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;
}

Approachable

For the record, because I had to search a little bit myself: The blitz manual 3.7.1 tells you how to promote user-defined types. Thanks for the hint. – Cannery

Answer 2

You should check this Wiki page: http://eigen.tuxfamily.org/dox-devel/TopicCustomizingEigen.html

Eigen requires some work to use other than primitive data types, but it's generally possible.

Cheep

Thanks for the hint. Had read the page and followed the hints. The point is, operator+ works just fine, but e.g. operator* is wrong, because meter * meter is not a meter. – Cannery

Answer 3

The difficulty of using the standard Eigen library plugin option, is that the existing operators +, -, *, etc, needs to be replaced for Boost Units quantities to be used.

For example, for a Boost units custom type to work with the * multiply operator, for an arbitrary CUSTOM_TYPE, it needs to look like this:

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>( ... );
}

Notice how the return type is not the same as the input types. Here you use the template helper multiply_typeof_helper to create the return type. This is because multiplying meters with seconds will not give you a quantity of either unit. However, the default Eigen * operator will return the same "type" as that of the inputs - this is the problem.

The other option is to embed the Eigen matrix inside the quantity, rather than embedding the quantity inside the matrix.

Chorea

Answer 4

The simple case of a vector / matrix with uniform units could partially be solved by adjusting the value type of an Eigen matrix, but this needs considerable user intervention in the form of specializing traits and even this is not sufficient for all operations.

However, if you want to solve the general case of matrices containing non-uniform units (i. e. each entry can have affigeren unit) then this will not work by customizing the value type of the matrix.

You need an approach like the one presented here which layers on top of a linear algebra library and provides the unit annotations: https://m.youtube.com/watch?v=4LmMwhM8ODI

Nano