Re: [eigen] Eigen and rigid body simulation

[Benoit Jacob <jacob.benoit.1@xxxxxxxxx>]

Right, your idea allows for much tighter integration with the rest of Eigen..
One of the main advantages of Maxime's proposal was that one could
make it a new module of its own; now with your idea, the whole
Geometry module depends on Lie stuff. That's not a bad thing: once we
agree that something is useful, better go 100% in that direction. Also
the basic Lie stuff doesn't add much code in itself. If there are more
specialized Lie things coming, like tangent spaces and whatnot, that
can still be content for a new "Lie" module.

Benoit

2009/12/2 Gael Guennebaud <gael.guennebaud@xxxxxxxxx>:
> 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;
>  };
>
> 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.
>
> 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 ?
>
> Also I don't understand why do you use functors for the log/exp functions
> and not static functions ?
>
> 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
>>
>>
>
>