Re: [eigen] Eigen and rigid body simulation

[Maxime Tournier <maxime.tournier@xxxxxxxxxxxx>]

On Wednesday 02 December 2009 13:52:17 Gael Guennebaud wrote:
> Yes that's what I'm thinking too, but now I'd like to hear Maxime's
>  opinion, in particular with such a design do you think it is still (I'm
>  quoting you):
> 
> "straightforward to implement the
> tangent bundle, then differentiable functions, and then automatically
> compute differentials for complicated function compositions using
> expression templates."
> 
> Since you have already done all of that in your framework, you're the best
> to tell us if the proposed design is not "wrong" :)
> 
> gael.
> 

I presume that the decorator approach could work for this too. 

On Lie groups the tangent bundle is trivial, that is isomorphic to a pair of 
one group element and one Lie algebra element (those are what you have to 
store internally).

basically you could implement it like: 

template<class G>
struct T {
  
  // the point of the base
  const G& at() const;

  // canonical isomorphisms: body and spatial velocities
  lie<G>::algebra body() const;
  lie<G>::algebra spatial() const;

  // named constructors
  static T body(const G& g, const lie<G>::algebra& v); 
  static T spatial(const lie<G>::algebra& v, const G& g);

 // implementation: one could only store the base point and e.g. the
 // body velocity, then the spatial velocity is obtained using the
 // adjoint map

 };

So yes, this seems compatible with the decorator approach.

To sum things up :

Decorator:
+ more readable
= perhaps more eigen-ish ? :)
- needs inheritance on types

Concept maps:
+ existing types untouched
- less readable

Maybe I can start with implementing a small, clean, documented, version of the 
code in unsupported/, so that others can start to play with, and based on the 
feedback we'll see what works in practice and what does not.


best,



max.


> On Wed, Dec 2, 2009 at 1:17 PM, Benoit Jacob 
<jacob.benoit.1@xxxxxxxxx>wrote:
> > 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
> 

-- 
Maxime Tournier

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