Go look up on
Wikipedia what it is. I am going to talk about how I am having to use it.
I was doing fine with the Rocketometer analysis code in C++, using the
NewMat library to handle matrices, with a Quaternion layer on top that I
wrote myself. After five days of doing battle with various things, I
finally got something that worked, but I was curious if this was the
"best" way to do it. The C++ Standard Template Library didn't have
anything directly related to matrices. The Boost library had a section
called uBLAS, but the documentation for it kind of de-recommended
itself. It suggested several alternatives, and the one that looked best
is called Eigen.
Eigen is interesting in that it is entirely header files, containing
almost all of its code in C++ templates. Templates are cool, mostly
because when they are instantiated, the compiler gets to see the code in
the context that it is used, and gets to inline and optimize it there.
Specifically, Eigen contains a dynamic-sized matrix, but also is a
template for fixed-sized vectors and matrices. I want to use these as
much as possible because all vector sizes used in Rocketometer data
analysis are known at compile-time, so the compiler can unroll loops and
so on to best optimize the code.
However, templates do not mix with virtual methods, so I had to figure
out how to make that work, since I used virtual methods to implement the
physics model. I had code that looks like this with NewMat:
class SimModel {
protected:
/** Physics function. Calculates the derivative of the state with respect to time
virtual ColumnVector fd_only(double t, const ColumnVector& x)=0;
/** Measurement function. Calculates the measurement from the given state and time */
virtual ColumnVector g_only (double t, const ColumnVector& x, int k)=0;
/** Physics function with process noise. Uses fd virtual function to calculate physics, then adds process noise.*/
ColumnVector fd(double t, const ColumnVector& x, const ColumnVector* v);
public:
/** Measurement function with measurement noise. Uses g virtual
function to calculate measurement, then adds measurement noise. */
ColumnVector g (double t, const ColumnVector& x, int k, const
ColumnVector* w) {ColumnVector result=g_only
(t,x,k);if(w)result+=*w;return result;};
};
But I wanted to adapt that to use Eigen, specifically with the
fixed-length vectors, since the size of the state vector is determined
by the problem and known at compile time. That means that ColumnVector
has to go, to be replaced by Matrix<double,n,1> where n is a
template parameter determining the size of the state vector. But what
about the measurement? The purpose of the k parameter to g_only is to
tell which of several kinds of measurements to use. For instance, in the
Rocketometer problem, we have a measurement vector coming from the
inertial and magnetic sensors, treated as a single 9-element vector. We
also have measurements coming from radar or equivalent, treated as a
3-element vector. So, we need a template function g_only, which
generates either a 9-element vector or a 3-element vector.
You can't do that and have it be virtual, too.
Basically, virtual functions are a runtime binding issue, while
templates are a compile-time binding. So, I can't have a virtual g_only
function, callable by the base class g function.
Enter the Curiously Repeating Template Pattern (CRTP). As it happens,
this is something that I read about just a few days ago, just reading up
on the C++ language in general. For us, the pattern goes something like
this:
template<int n, class Derived> class SimModel {
public:
template<int m> Matrix<double,m,1> g (double t, const Matrix<double,n,1>& x, int k, const Matrix<double,m,1>* w) {
Matrix<double,m,1> result=static_cast<Derived*>(this)->template g_only<m>(t,x,k);
if(w)result+=*w;
return result;
};
};
Note that g_only
isn't even defined in this class template, only
used. In fact, one of the weaknesses of CRTP is that it implies
definitions without expressing them, so it is hard to document. Also note that extra template keyword there after the arrow operator. See
here for details.
The
derived class then looks like this:
template<int n> class RocketometerModel: public SimModel<n,class RocketometerModel<n>> {
template<int m> Matrix<double,m,1> g_only(double t, const Matrix<double,n,1>& x, int k);
};
So what happens is that the compiler
- Parses the template for SimModel, but doesn't compile it, because
it's a template, not actual code yet. Therefore it doesn't matter that
g_only is undefined yet.
- Parses the template for RocketometerModel, and again doesn't compile it.
- Parses the main code, compiling as it goes along until it hits
RocketometerModel. It instantiates and compiles RocketometerModel, in
the process instantiating and compiling SimModel.
- When SimModel is being instantiated and compiled, it has a call to
RocketometerModel g_only, but this is ok, since that is available
already, since step 2 has already happened.
Now the derived class might not be a template, it might be an actual
class. In this case, the base class is instantiated and compiled with
the derived class. In either case, everything is available before it is
used, even though the code might look otherwise.
Now this part I will write in bold, so that Google can see it. The
other curiously repeating template pattern is having to use the word .template when using (not defining) a template member function.
This solves the error invalid operands of types ‘<unresolved overloaded function type>’ and ‘int’ to binary ‘operator<’ .
It appears that there is a weakness in the C++ operator precedence table, such that in certain cases it can't distinguish the opening angle bracket of a template parameter set from a normal compare-for-less-than. In order to disambiguate, you throw in the word template after the dot (it will also work after a pointer arrow -> if you are using one of those). I don't understand it completely myself, but Eigen uses this itself, which is how I found out about how it works in the first place.
I am gradually coming to the conclusion that Java did it right, with generics which are real compiled code. Java generics are enabled by the "everything is an object" model in which all objects descend from a common class. I am also beginning to think that Java did it right with operator overloading. Operator overloading, even when fully appropriate, like defining * to mean matrix or quaternion multiplication, is fun to
use but a nightmare to
implement. And, if it is implemented wrong, it might be left to the user to find out when he tries to do something the implementer did not foresee.
All in all, I give Eigen 4/5, especially recommended for new projects, but not for converting old projects. The biggest advantage is speed. What took IDL over an hour, took Java and C++ with NewMat about 4 minutes, but takes Eigen only 20 seconds. Also, templated matrix and vector sizes are nice, because they resolve matrix size mismatches at compile-time. Finally, zero-based component indexing is what I expect, and the reason I don't suggest converting old projects from NewMat. Also be aware that the Eigen quaternion library uses the convention \(\vec{v}_{to}=\mathbf{p}\vec{v}_{from}\mathbf{p}'\), which is fine and internally consistent, but not consistent with the math I had seen for quaternion derivatives. As a consequence, my code is liberally festooned with q.conjugate() and in some places q.conjugate().conjugate(). It's
almost a case of two wrongs making a right, but not quite.
