Implement differently sized-pipelines in Expression-inl.h

[Frank Dellaert]

Currently, only the m=2 pipeline is implemented. It's a drag to add new pipelines, however, because there is a lot of copy/paste. At one point, @hannessommer , you sent email about a scheme to make this be more or less automatic, but it was still somewhat elaborate. Might you have time to look at Expression again and see what your thoughts are now?

This is a critical enhancement as it has a large effect on performance.

[Hannes Sommer]

Ah, yes. I can certainly convert the ideas from back then into a revision of the current code. Probably this WE. There is a full automatic solution: one would only need to provide a compile time integer to specify the limit MAX_STATIC_ROWS above which dynamic row numbers would be used. But that would involve having some classes of MAX_STATIC_ROWS as inheritance depth. If that is no problem I could do that - for the compilers that is totally fine. Otherwise one has to write all the MAX_STATIC_ROWS many virtual function declarations and once all its implementations by a class template (only one line each). That is then the least structural duplication with inheritance depth constraints.

[Frank Dellaert]

Awesome. Looking at it again, the whole things starts in

    virtual void XXXNode::Record::startReverseAD(JacobianMap& jacobians) const {
      Base::Record::startReverseAD(jacobians);
      Select<traits::dimension<T>::value, A>::reverseAD(This::trace, This::dTdA, jacobians);
    }

but, I'm now wondering whether we need a limit at all? I think maybe other parts of the pipeline already make it impossible to have a dynamically-sized return type. Hence, why bother with the dynamic pipeline at all? It will just clutter the code...

[Hannes Sommer]

The problem back then was that you would need to specify virtual functions for each number of rows. And you can not have them added later on demand. Basically the problem of not having templated virtual functions ;). So one has to specify which static row sizes are supported. To catch any size of matrix one needs another virtual function taking an dynamically sized matrix, right?

[Frank Dellaert]

I remember it being a brilliant hack :-) to what I then thought was indeed a problem like that, but looking at the code now I have no clue why

    virtual void Record::startReverseAD(JacobianMap& jacobians) const {
      Base::Record::startReverseAD(jacobians);
      Select<traits::dimension<T>::value, A>::reverseAD(This::trace, This::dTdA,
          jacobians);
    }

can't simply call a templated version of

  typedef Eigen::Matrix<double, 2, Dim> Jacobian2T;
  void ExecutionTrace::reverseAD2(const Jacobian2T& dTdA, JacobianMap& jacobians) const {
    if (kind == Leaf)
      handleLeafCase(dTdA, jacobians, content.key);
    else if (kind == Function)
      content.ptr->reverseAD2(dTdA, jacobians);
  }

Record knows both ROWS== traits::dimension<T>::value and COLS== traits::dimension<T>::value, but they must be different T's :-) This is a complex thing to keep on top of, fading fast :-)

[Hannes Sommer]

Yes and no, at that particular place you are right and the two functions :

  // Either add to Jacobians (Leaf) or propagate (Function)
  void reverseAD(const Matrix& dTdA, JacobianMap& jacobians) const {
    if (kind == Leaf)
      handleLeafCase(dTdA, jacobians, content.key);
    else if (kind == Function)
      content.ptr->reverseAD(dTdA, jacobians);
  }
  // Either add to Jacobians (Leaf) or propagate (Function)
  typedef Eigen::Matrix<double, 2, Dim> Jacobian2T;
  void reverseAD2(const Jacobian2T& dTdA, JacobianMap& jacobians) const {
    if (kind == Leaf)
      handleLeafCase(dTdA, jacobians, content.key);
    else if (kind == Function)
      content.ptr->reverseAD2(dTdA, jacobians);
  }

In ExecutionTrace can be just replaced by a templated function and treat every matrix statically, but internally this function has to do call reverseAD for the content.ptr. This needs would be a virtual function call. Hence, the function called then can't be a templated function. So it needs to pick a virutal function based on the dimension (either by some specialized template, as it is now, or by overload resolution). And you cannot completely get rid of this one step of virtual function call unless you do fully typed expressions (the whole expression tree is part of the type). There are also all kind of compromises between the two. (clusters of complex expression types linked through virtual function calls)

[Frank Dellaert]

Right! OK...

[Hannes Sommer]

First proposal in pull request #38

[Frank Dellaert]

Hannes nailed it :-)

[Andrew Melim]

I'm reopening this issue since the current implementation does not compile on windows. I'm taking a look at it now.

[Frank Dellaert]

What error are you getting and where?

[Andrew Melim]

Error   22  error C2782: 'void gtsam::GenerateFunctionalNode<U1,U2,U3>::Record::reverseAD(const Eigen::Matrix<double,Rows,Cols,0|_Rows==1&&_Cols!=1?RowMajor:_Cols==1&&_Rows!=1?0:0,_Rows,_Cols> &,gtsam::JacobianMap &) const' : template parameter '_Rows' is ambiguous   C:\Users\water_000\borg\gtsam\gtsam_unstable\nonlinear\CallRecord.h 106 1   Pose2SLAMExampleExpressions

[Andrew Melim]

I'm looking through the code and I'm not completely clear on what is going on :/ Not comfortable hacking away at getting it to compile. Any thoughts?

[Frank Dellaert]

Even I am lost in this code now, since @hannessommer went to town on this. @hannessommer , I know this has probably nothing to do with the recursive calls, but (a) do you have a clue on how to fix on windows? (b) would you be able to do a pull request where you manually unroll, as ou suggested might be more readable?

[Hannes Sommer]

I will have a look tomorrow morning. Looks like a buggy compiler :) - which one?

[Andrew Melim]

VS2013-v4

[Hannes Sommer]

Looking at this one line of error. This is complete wired. According to it a template parameter ambiguity is supposed to happen when the body of this (non-templated) function:

  virtual void _reverseAD(
      const Eigen::Matrix<double, MaxSupportedStaticRows, Cols> & dFdT,
      JacobianMap& jacobians) const {
    static_cast<const Derived *>(this)->reverseAD(dFdT, jacobians);
  }

Calls this templated function:

    /// Given df/dT, multiply in dT/dA and continue reverse AD process
    template<int Rows, int Cols>
    void reverseAD(const Eigen::Matrix<double, Rows, Cols> & dFdT,
        JacobianMap& jacobians) const {
      Base::Record::reverseAD(dFdT, jacobians);
      This::trace.reverseAD(dFdT * This::dTdA, jacobians);
    }

This is nonsense to me :). @amelim, can you provide more of the error message and maybe the commit id you are using? And unrolling would not help here. But still I can do this. Maybe at the weekend.

[Frank Dellaert]

Maybe this is another case of the options: I have seen Windows complain about that earlier. @amelim , also quickly try changing argument to `Matrix<double, Rows, Cols, Options, _Rows, _Cols>

[Frank Dellaert]

And, @hannessommer , unrolling would be great. I think we went a bridge too far - even thinking of undoing the first recursive construction. It is probably opaque to anyone but us, and does not save a lot of space.

[Hannes Sommer]

Generalizing the signature of the template function could be a good idea. I would even go for:

  template<typename Derived>
  inline void reverseAD(const Eigen::MatrixBase<Derived> & dFdT,
      JacobianMap& jacobians) const

This is even simpler for the compiler.

[Hannes Sommer]

pull request #56 would be the unrolling. But I don't believe this step will help the windows compiler.

[Frank Dellaert]

@amelim please try this again

[Andrew Melim]

Took a look at this again today and the issue is still not resolved. I'm working through the suggestions above.

[Andrew Melim]

The issue is now isolated to the following

Error   3   error C2784: 'void gtsam::GenerateFunctionalNode<U1,U2,U3>::Record::reverseAD4(const Eigen::Matrix<double,Rows,Cols,0|_Rows==1&&_Cols!=1?RowMajor:_Cols==1&&_Rows!=1?0:0,_Rows,_Cols> &,gtsam::JacobianMap &) const' : could not deduce template argument for 'const Eigen::Matrix<double,Rows,Cols,0|_Rows==1&&_Cols!=1?RowMajor:_Cols==1&&_Rows!=1?0:0,_Rows,_Cols> &' from 'const gtsam::Matrix' C:\Users\water_000\borg\gtsam\gtsam_unstable\nonlinear\CallRecord.h 165 1   Pose2SLAMExampleExpressions

[Frank Dellaert]

Try replacing with

    template<typename SomeMatrix>
    void reverseAD4(const SomeMatrix & dFdT, JacobianMap& jacobians) const {

[Andrew Melim]

That fixed it!

[Frank Dellaert]

Yay !