How to integrate a library that uses expression templates?

Question

I would like to use the Eigen matrix library as the linear algebra engine in my program. Eigen uses expression templates to implement lazy evaluation and to simplify loops and calculations.

For example:

#include<Eigen/Core>

int main()
{
  int size = 40;
  // VectorXf is a vector of floats, with dynamic size.
  Eigen::VectorXf u(size), v(size), w(size), z(size);
  u = 2*v + w + 0.2*z;
}

Since Eigen uses expression templates, code like

u = 2*v + w + 0.2*z;

In the above mentioned sample reduces to a single loop of length 10 (not 40, the floats are put into regiser by chunks of 4) without creating a temporary. How cool is that?

But if I integrate the library like this:

class UsingEigen
{
    public:  
        UsingEigen(const Eigen::VectorXf& data):
            data_(data)
        {}

        UsingEigen operator + (const UsingEigen& adee)const
        {
            return UsingEigen(data_ + adee.data_);
        }

        ...
    private:
        Eigen::VectorXf data_;
}

Then the expressions like:

UsingEigen a, b, c, d;
a = b + c + d;

cannot take advantage of the way Eigen is implemented. And this is not the last of it. There are many other examples, where expression templates are used in Eigen.

The easy solution would be not to define the operators by myself, make data_ public and just write expressions like:

UsingEigen a, b, c, d;
a.data_ = b.data_ + c.data_ + d.data_;

This breaks encapsulation, but it preserves the efficiency of Eigen.

Other way could be to make my own operators, but let them return expression templates. But since I am a beginner in C++, I do not know if this is the right way to go.

I am sorry if the question is too general in nature. I am a beginner and have noone to ask. Up until now I was using std::vector<float> everywhere, but now I need to use matrices also. To switch from std::vector<float> to Eigen in my entire project is a big step and I am afraid of making a wrong call right in the start. Any advice is welcomed!

Vanderhoek

Interesting problem. The first thing, though, is: why do you want to encapsulate the Eigen library vectors in this way? What behaviour do your classes add? – Ormuz

My classes do not add any functionality to the Eigen library classes themselves, but use them. For example one of my core classes stores two vectors. One as an input to a certain mathematical computation and the other as output. These objects need to interoperate in a similar fashion as I mentioned above. When you add two such objects, the inputs should be added. – Vanderhoek

I don't think this is possible without reproducting substantial part of expression template framework yourself. For instance, (a+b)*c will be something like ExprCwiseAdd*UsingEigen (the name is made up, don't recall it any longer), and there will have to be ExprCwiseAdd*UsingEigen defined somewhere, but also ExprCwiseAdd*ExprCWiseAdd and so on. In short, the addition will not have UsingEigen as return type. (You might have a look at boost::proto which is a framework for expression templates). Good luck. – Lilliamlillian

Also, related to this topics are two tenets I learned from books: 1. there is rarely any good reason to make data public. 2. When using external libraries, rest of the project should be shielded from changes in them using an interface class. – Vanderhoek

@Martin In that case I would try to model all parts of the workflow that do perform calculations purely in terms of Eigen classes, and only use your wrappers as the respective input and output. Intermediate steps then use the actual Eigen data. And you don’t need to expose variables, you can use functions that return (for instance) references to the classes. True, this still exposes data. I feel your pain. – Ormuz

If you ever decide to use expression templates in your own code, do yourself a favor and take a look at Boost.Proto. – Cordes

Answer 1

Why would exposing data_ break encapsulation? Encapsulation means hiding the implementation details and only exposing the interface. If your wrapper class UsingEigen does not add any behavior or state to the native Eigen library, the interface does not change. In this case, you should drop this wrapper altogether and write your program using the Eigen data structures.

Exposing a matrix or a vector is not breaking encapsulation: only exposing the implementation of the matrix or vector would do that. The Eigen library exposes the arithmetic operators but not their implementation.

With expression template libraries, the most common way for users to extend the library functionality is by adding behavior, not adding by adding state. And for adding behavior you do not need to write wrapper classes: you can also add non-member functions that are implemented in terms of the Eigen class member functions. See this column "How Non-Member Functions Improve Encapsulation" by Scott Meyers.

As for your concern that the transformation of your current program to a version that explicitly uses the Eigen functionality: you can perform the change step-by-step, changing small parts of your program each time, making sure your unit tests (you do have unit tests, don't you?) do not break as you go along.

Sather

Exposing data members breaks encapsulation. But adding the expression templates overloads for the Eigen classes would break it in the same way. – Ormuz

@KonradRudolph It doesn't break encapsulation if UsingEigen has the exact same interface and with the exact same implementation (i.e. a literal wrapper class with pure forwarding and no logging/checking etc. added). You would be right if UsingEigen would define an extra layer of abstraction, but it doesn't appear to be the case (and it should not be called UsingEigen then at all because that exposes the implementation!) – Sather

Answer 2

In my opinion, this looks more of a object oriented design problem rather than a library usage problem. Whatever you read from the books are the right recommendations. i.e, do not expose member variables and shield the upper layers from the nuances of the 3rd party layer usage.

What you could look forward is right abstractions of mathematical functions that can be implemented using this library internally. i.e, you could expose a library of your own with high level functions than elementary vector and matrix operations. In this way you can utilize the peculiarities of the interactions among the library objects and at the same time you don't have to expose your member variables to upper layers.

For e.g you could abstract away my higher level APIs like computing the distance from a point to a plane, distance between two planes, computing the new coordinates of a point w.r.t another coordinate system using the transformation matrices etc. To implement these methods internally you can utilize the library objects. You can restrict to not to have any of the library classes used in the API signatures to avoid dependency for the upper layers on this library.

Upper layers of your program should be higher in the level of abstraction and need not bother about the elementary implementation details such as how the calculation of the distance from a point to the plane is implemented etc. Also, they even need not know if this lower layer is implemented using this library or something else. They would just use the interfaces of your library.

Gaspar

Answer 3

Set up a class template to hold general Eigen expressions and make UsingEigen a special instance of it:

template<typename expr_t>
class UsingEigenExpr
{
    UsingEigen(expr_t const& expr) : expr(expr) {}
    expr_t expr;
    operator UsingEigenExpr<Eigen::VectorXf>() const
    {
        return {expr};
    }
};
using UsingEigen = UsingEigenExpr<Eigen::VectorXf>;

Then overload any required function, e.g. as

template<typename expr1_t, typename expr2_t, typename function_t>
auto binary_op(UsingEigenExpr<expr1_t> const& x, UsingEigenExpr<expr2_t> const& y, function_t function)
{
    return UsingEigenExpr<decltype(function(std::declval<expr1_t>(),std::declval<expr2_t>()))>(function(x.expr,y.expr));
}

template<typename expr1_t, typename expr2_t>
auto operator+(UsingEigenExpr<expr1_t> const& x, UsingEigenExpr<expr2_t> const& y)
{
    return binary_op(x,y,[](auto const& x, auto const& y) {return x+y;});
}

and so on for other binary operators like operator-, for unary operators, and more generally for all the other stuff you want to use. Further, you could add some other member functions to UsingEigenExpr, e.g. size(), norm(), etc.

Use it as

UsingEigen b, c, d;
auto a = b + c + d;

to store the expression, or

UsingEigen b, c, d;
UsingEigen a = b + c + d;

to directly evaluate it.

Although this approach works, in the end you find yourself duplicating all the required functionality, so use it carefully.

Quillet

Answer 4

I don't understand all your question I will try to answer you most of them. In this sentence:

 UsingEigen operator + (const UsingEigen& adee)const
    {
        return UsingEigen(data_ + adee.data_);
    }

You have an overload operator ( sorry I don't know if this is the correct way to write in English), for this reason you can write:

a = b + c + d;

instead of:

a.data_ = b.data_ + c.data_ + d.data_;

You won't have any problem, cost of your program will be the same. In addition you will have encapsulation and efficiency.

On the other way if you want define your own operator you can do it like the template do it. You can find information on the web searching "overload operator" but is similar to this:

UsingEigen operator + (const UsingEigen& adee)const
    {
        return UsingEigen(data_ + adee.data_);
    }

Instead of "+" you can put the operator and do the operations you need.

If you want to create a matrix it is simple. You only need to create a array of array or vector of vector.

I think is something like this:

std::vector<vector<float>>

I am not sure but it is easy, on the other hand you can use a simple matrix on this way:

float YourMatrix [size][size];

I hope it could help you. I don't understand all your question if you need something more add me on google+ and I will try to help you.

Sorry for my English, I hope you can understand all and it helps you.

Guncotton

Unfortunately, you are wrong. You need to read up on what expression templates are and how they influence this code. Your overloaded operators won’t take advantage of expression templates, exactly as the OP has described. Furthermore, the Eigen library provides way more semantics than the naive nested-array approach you’ve shown, which is completely inadequate for OP. – Ormuz

Thank you Manuel, If I use your approach, the program will compile and work correctly. But expressions like a = b + c + d; will result in something like this: first b + c will be computed and stored in a temporary bctemp. Then bctemp will be added to d and create bcdtemp. At last, the bcdtemp will be assigned to a. That is 3 loops and 2 temporaries. My goal is to avoid this. – Vanderhoek

@Manuel If you are interested in what I mean by avoiding the temporaries and loops, you can find better explanation here eigen.tuxfamily.org/dox/TopicInsideEigenExample.html – Vanderhoek

thank you!I will read it. I am sorry for not be able to help you. :( I hope somebody help you as soon as possible. – Guncotton