Re: [eigen] Eigen and rigid body simulation

[Maxime Tournier <maxime.tournier@xxxxxxxxxxxx>]

Wow, lots of questions to answer :)

On Wednesday 02 December 2009 11:31:03 Gael Guennebaud wrote:
> I also like Maxime's idea, but I'm wondering whether the Lie<G> group class
> should not be a decorator to G rather than a traits class. The main raisons
> for that are 1) to make it less tedious to use, 2) the current design of
>  the geometry module already goes to that direction since all our
>  transformation classes already implement inverse, identity and composition
>  via operator* (e.g., Translation::operator* actually performs a +).
> 
> So basically, a Lie object should implement the following interface:
> 
> template<class G>
>   struct Lie {
>   typedef some_type Algebra;
>   typedef some_type CoAlgebra;
> 
>   static G Identity();
>   G inverse() const;
>   G operator*(const G&) const;
> 
>   G exp() const;
>   G log() const;
>  };
> 

yep, you could do it like this, but then you have to take care about storage 
in each specialization, and struggle with name clashes introduced by 
operators, exp/log, and so on.

Also, and this is really important, one *has* to derive from Lie<G> in order 
to use the Lie structure. This is a strong constraint, 'cause now it's 
impossible to say that a pair of group elements is actually an element of some 
product group, unless you find a way to make std::pair<A, B> derive from Lie< 
std::pair< A, B > > which is of course impossible.

This is the very motivation of using "trait classes": it leaves existing types 
unchanged. This will make it much easier for people to benefit from the 
framework since all they have to do is to provide a concept_map (in c++0x's 
wording) for the concept Lie, and everything should work. They wont have to 
mess with their oh-so-shiny class hierarchy to make it work. It'll just stay 
shiny.

Of course, the burden is somehow shifted towards library developpers (=us) but 
hey, this is gonna be a template nightmare anyways :-)

> and then, e.g., Translation<S,D> could inherit Lie<Translation<S,D> >,
> NonUniformScaling<S,D> would inherit Lie<NonUniformScaling<S,D> >,
> Quaternion<S> would inherit Lie<Quaternion<S>>, etc.
> 

> Each of these specializations of Lie<> would also provide the storage and
> respective ctors. This way one can still directly instanciate a
> Lie<Translation<S,D> > if he really wants a "pure" Lie group object, and
> more importantly, write generic algorithms taking only objects implementing
> the Lie group concept:
> 
> template<G>
> Lie<G> fast_exp(const Lie<G>& x, int n)
> {
>   if( n == 0 )
>     return Lie<G>::Identity();
> 
>   if( n % 2 )
>     return x * fast_exp(g*g, n/2);
>   else
>     return fast_exp( g*g, n/2);
> }
> 
> that is much more readable.
> 
totally agreed. 

and let's not even mention those spurious "typename" compiler hints. 

> However, the exp/log functions could now become ambiguous, especially for
> Quaternion. So if you think the above proposal makes sens, we could use
> lieExp/lieLog to make it more explicit they belong to the Lie algebra.
> 
> 
> Now I have to acknowledge that your design has so advantages too. In
> particular it really well separate the Lie algebra, and make it clear when
> you are using operations belonging to the Lie algebra. In that sense it
> looks a bit cleaner. Are there any other advantages that I'm missing ?
> 
in a nutshell: it doesnt touch existing types.

> Also I don't understand why do you use functors for the log/exp functions
> and not static functions ?
> 

i could have done it with static functions as well (it used to be the case 
actually), but i also happen to derive these functions which is far easier 
with functors: the functors have a .diff method from the domain tangent bundle 
to the range tangent bundle. 

Actually i found functors a lot simpler to deal with than plain functions, 
once you get used to the (small) extra coding.

Benoit, as for the relationships between this and what Mathieu is doing, from 
my point of view it is the very same thing in the end (after all, the goal is 
to simulate rigid bodies in both cases right ?), using two different 
approaches.

Mathieu's one is perhaps more user-friendly while mine seems more developper 
friendly (at least to me).

For example I don't see why one should introduce (and maintain) a Twist class 
since this essentially boils down to a 3-vector, same with Torque (3-covector) 
or Wrench (6-covector). All you basically need is an adjoint/coadjoint 
operator on these spaces.



max.


> gael.
> 
> On Tue, Dec 1, 2009 at 11:24 PM, Benoit Jacob 
<jacob.benoit.1@xxxxxxxxx>wrote:
> > 2009/11/30 Maxime Tournier <maxime.tournier@xxxxxxxxxxxx>:
> > > Hi everyone,
> > >
> > > As Mathieu said, I have some generic code that allows one to abstract
> > > the Lie group structure of a type, so that one can write any Lie group
> > > algorithm no matter the group. I told Gael about it few days ago, and
> > > since it could be of interest to others, here's an overview of what I
> > > did.
> > >
> > > The approach somehow follows c++0x's concepts proposal, in that the
> > > Lie group structure for a type T is encoded in a separate template
> > > class, Lie<T>, containing all the types and operations that define the
> > > Lie group structure.
> > >
> > > One can describe the Lie group structure for a type T by specializing
> > > the template Lie for this type, giving the associated types and
> > > implementing required operations.  Every specialization should conform
> > > to the same "interface".
> > >
> > > Here are two basic examples with only the algebraic structure:
> > >
> > > // here is S^3 with quaternion multiplication
> > > template<class Real>
> > > struct Lie< Quaternion<Real> > {
> > >
> > >  static Quaternion<Real> id() { return Quaternion<Real>::Identity(); }
> > >  static Quaternion<Real> inv(const Quaternion<Real>& q) { return
> >
> > q.inverse();
> >
> > > }
> > >  static Quaternion<Real> comp(const Quaternion<Real>& q1,
> > >                               const Quaternion<Real>& q2) {
> > >    return Quaternion<Real>::Identity();
> > >  }
> > >
> > > };
> > >
> > > // here is R^3 with addition
> > > template<>
> > > struct Lie< Vector3d > {
> > >
> > >  static Vector3d id() { return Vector3d::Zero(); }
> > >  static Vector3d inv(const Vector3d& v) {
> > >    return -v;
> > >  }
> > >  static Vector3d comp(const Vector3d& a,
> > >                       const Vector3d& b) {
> > >    return a + b;
> > >  }
> > >
> > > };
> > >
> > > You can then write an algorithm for any kind of group like this:
> > >
> > > template<class G>
> > > G fast_exponentiation(const G& g, int n) {
> > >  if( n == 0 ) {
> > >    return Lie<G>::id();
> > >  }
> > >
> > >  if( n % 2 ) {
> > >    return Lie<G>::comp(  fast_exponentiation( Lie<G>::comp(g, g), n/2),
> > > g
> >
> > );
> >
> > >  } else {
> > >    return fast_exponentiation( Lie<G>::comp(g, g), n/2) );
> > >  }
> > >
> > > }
> > >
> > > Now I know it's a bit tedious to write, but:
> > >
> > > - No modification to exisiting types is needed whatsoever. This is A
> > >  Good Thing. No curious inheritance pattern is introduced in
> > >  particular. Actually, underlying types are completely independant
> > >  from this new structure. This should make Eigen types developers
> > >  happy since we're not touching anything in existing Eigen code :-)
> > >
> > > - We can quite easily construct new groups from existing ones:
> > >
> > > template<class G1, class G2>
> > > struct Lie< std::pair<G1, G2> > {
> > >
> > >  // ...
> > >
> > > };
> > >
> > > ...which essentialy says that the product of two groups is a group for
> > > the product law. This is particularly useful when using std::tuple and
> > > variadic templates, since product-groups are automatically constructed
> > > at compilation.
> >
> > I understand, your design seems very good.
> >
> > > The complete Lie template specification in my code is the following:
> > >
> > > template<class G>
> > > struct Lie {
> > >  typedef some_type algebra;
> > >  typedef some_type coalgebra;
> > >
> > >  static G id();
> > >  static G inv(const G&);
> > >  static G comp(const G&, const G& );
> > >
> > >  // default-constructible
> > >  typedef some_functor exp;
> > >  typedef some_functor log;
> > >
> > >  // G-constructible
> > >  typedef some_functor ad;
> > >  typedef some_functor ad_T;
> > >
> > > };
> >
> > OK I understand. ad is inner automorphisms? and ad_T is the inverse?
> >
> > > The only problem with it is the identity for dynamic-sized types, for
> > > example VectorXd: what should be the size of the result ?
> >
> > Do you mean as a compile time value or as a runtime value?
> >
> > At compile time:
> >
> > Throughout Eigen we use a special value, Dynamic, to mean that. For
> > example, VectorXd is just:
> > matrix<double,Dynamic,1>.
> > Is that what you were looking for?
> >
> > At run time: in the dynamic-size case, you could require id() to take
> > a size argument. So id(void) would use a static assertion that the
> > size is fixed at compile time. See what is already done in
> > MatrixBase::Identity().
> >
> > > So maybe it
> > > should accept some optional contextual hint for the dimension.
> >
> > in Eigen, all types inheriting MatrixBase (hence VectorXd) have a
> > member enum SizeAtCompileTime telling the number of coefficients
> > (rows*cols)
> >
> > so:
> >   VectorXd::SizeAtCompileTime == Dynamic
> >
> > Other types could be harmonized if that is useful to you, e.g. we
> > could set Quaternion::SizeAtCompileTime==4 if that's useful.
> >
> > > Once this is done, it is quite straightforward to implement the
> > > tangent bundle, then differentiable functions, and then automatically
> > > compute differentials for complicated function compositions using
> > > expression templates. When dealing with complicated functions on
> > > complicated non-flat spaces, this is a real joy to use.
> > >
> > > I have the Lie group structure implemented for quaternions, Rohit
> > > Garg's rigid transforms, row/column vectors, as well as any std::tuple
> > > combination of these. Based on this I have code that can interpolate
> > > (linear/~spline), compute geometric mean, and do statistical analysis
> > > while preserving the Lie group (hence manifold) structure.
> > >
> > > So if you need to spline interpolate homographies
> > > in a non degenerate way for example, all you have to do is to
> > > implement the Lie group structure for homographies and the algorithm
> > > is already coded for you :-)
> > >
> > > So my question is: would anyone be interested in me porting all/parts
> > > of this stuff into Eigen ? If so where/how should I start ?
> >
> > It seems that this would be nice content for a new module, "Lie". Does
> > that sound good?
> >
> > You would start by forking the eigen repository at bitbucket, and
> > creating a new Lie module in unsupported/.
> >
> > If you prefer it's always possible to develop directly in our
> > repository, it's just that forking gives you that extra flexibility.
> >
> > I would also like to understand how this interplays with Mathieu's
> > proposal. If you create a new Lie module, should Mathieu's classes go
> > into it? Or should only Torque/Twist/Wrench go into it while "so(3)"
> > is considered so important for geometry that it goes into the Geometry
> > module? If Mathieu's stuff gets split between Geometry and Lie, does
> > that mean that Lie should #include Geometry? Is that OK? Is that
> > something you wanted to do anyway?
> >
> > Benoit
> 

-- 
Maxime Tournier

PhD Student - EVASION Team
Grenoble Universities / LJK / INRIA Rhône-Alpes
FRANCE
Tel: +33 4 76 61 54 45