Auto-generate Matlab bindings for our C++ libraries

[rdeits]

Long-term, it would be nice not to have to write every mexfunction by hand.

SWIG is a mature tool designed to produce bindings for a C/C++ library in the higher-level language of your choice (python, java, perl, R, etc.). I've been playing with it for a few days using Python, and it seems really nice (see: https://github.com/rdeits/swig-eigen-numpy ). It handles a lot of the things that we've spent a whole lot of programmer time on, like wrapping pointers to C++ objects and converting data types to and from their C++ equivalents. You can create custom type conversions (like automatically converting numpy arrays in python to Eigen Matrices in C++), and that seems to work really well.

Unfortunately, it doesn't support MATLAB. There's an experimental version of swig which claims to have basic MATLAB support: http://stackoverflow.com/a/26853963 , but I haven't tried it yet.

(SWIG does support Octave, which looks to have something kind of like a mexfunction interface: http://www.swig.org/Doc3.0/Octave.html#Octave so that might be another possible route to take. But given that the CasADi people have been working on swig-matlab for over a year now, I imagine that they've decided the Octave route isn't a good idea).

I still think that something like this could be a really big improvement for our lives. Things like DrakeMexPointer work now but have been a huge timesink for us in the past. Maybe somebody in the group wants to finish up the SWIG-MATLAB interface ;-)

[RussTedrake]

i took a stab at what we might want this to look like (by generating an example manually).

First observation: everything is cleaner if the matlab class is a handle class. changing the DynamicalSystems hierarchy to handle classes is a big interface change. but it might be the right thing to do for the long run.

Therefore, I started with CoordinateFrame. take a look at this:
https://github.com/RussTedrake/drake/tree/coordinate_frame_cpp/drake/systems/frames/%40CoordinateFrame
I'm pretty happy with it. It means generating a LOT of mex wrapper files for more complicated classes. But they're all small. If they get generated automagically by swig, then it seems reasonable. And the interface it provides is pretty clean.

Notice that I converted part of the functionality of CoordinateFrames to c++, but not all. Partly because I couldn't, but also because I wanted to see how it might look when the matlab classes implement extra stuff on top of the c++ implementation. I think it's pretty clean.

I'm running unit tests locally now, but might actually push this. I'm curious what people think.

[patmarion]

Is it better to let matlab dispatch to your mexfunctions (one mexfunction per c++ method, as demonstrated here) versus doing the dispatch yourself (one mexfunction per class, first arg is a method name as a string)?

[psiorx]

I'd say it's better to have matlab dispatch to the mexfunctions. We tried it the other way a while back with RigidBodyManipulatormex and it turned to to be both slower and more unwieldy to work with.

[RussTedrake]

I've changed my mind about this one. The ability to have c++ bindings appear in the matlab namespace by simply installing to the class directory is significantly easier and cleaner than having to encode on the matlab side and decode on the mex entry side.

On Aug 19, 2015, at 11:43 AM, John Carter notifications@github.com wrote:

I'd say it's better to have matlab dispatch to the mexfunctions. We tried it the other way a while back with RigidBodyManipulatormex and it turned to to be both slower and more unwieldy to work with.

—
Reply to this email directly or view it on GitHub.

[RussTedrake]

I love the idea that we could use swig to generate bindings for many languages -- not just matlab. But one concern is that it is gpl 3. We should be careful to make sure we use it in a way that it does not infect drake. That might be slightly harder if we have to write our own language bindings...

[rdeits]

Looks like the license situation is fine: "When SWIG is used as it is distributed by the SWIG developers, its output is not governed by SWIG's license (including the GPL)." More detail here: http://www.swig.org/legal.html

[RussTedrake]

Regarding the speed of the many entry point dispatch - i believe that the first call to a method is relatively slower (matlab will dynamicaly load from disk) but that repeated calls should not be any slower -- assuming things aren't getting paged out.

@psiorx - is that consistent with what you remember?

[RussTedrake]

Right, but what if we have to write the matlab bindings?

On Aug 19, 2015, at 12:17 PM, Robin Deits notifications@github.com wrote:

Looks like the license situation is fine: "When SWIG is used as it is distributed by the SWIG developers, its output is not governed by SWIG's license (including the GPL)." More detail here: http://www.swig.org/legal.html

—
Reply to this email directly or view it on GitHub.

[rdeits]

As I see it we'd probably have our own fork of swig based on their matlab branch, which would be GPLed, but keep it separate from Drake. In that case we can keep Drake under whatever license we want, since it would be just another user of (in this case, our version of) swig.

They also specifically mention that some parts of the swig codebase get copied directly into output files, and those parts are licensed to be redistributed without restriction. So even if some swig code ends up in the drake output, we can still use the license of our choice.

[psiorx]

@RussTedrake
Yes, the first call will be slower. Once the mex file is resident in memory, subsequent calls will be normal speed.

[RussTedrake]

Summarizing our conversations in the lab:

• more mex entry points is better. @psiorx and I agree.
• we modify swig as necessary and keep it as an external. drake�s license will not be infected.

[rdeits]

I've investigated the state of the swig-matlab project, with mixed results. Wrapping C/C++ functions seems to work well, and swig even puts all the bindings in a nice Matlab package for you. It's also possible to get and set global variables inside the C++ file without writing your own getters and setters. Classes also work, and we can instantiate them, modify their properties (again, without writing our own getters and setters), and call methods.

But there are also some downsides. The most significant is that there seems to be a bug in the way their class deletion wrappers are written, which results in your objects never actually being destructed. I've got an open issue to investigate that: jaeandersson/swig#36 . Also, swig-matlab generates a single mex function which then dispatches to all of your C++ methods, which we already decided we weren't in favor of. They do, however, generate Matlab stubs for every function, so tab-completion does work fine.

I would guess that our DrakeMexPointer logic is probably more mature than their equivalent system.

[tkoolen]

Late to the party, but I came up with what I think is a really cool way to standardize mexification of functions. Basically, it would reduce the body of every mex function to two lines.

I did a quick mockup where I used this approach to mexify massMatrix and geometricJacobian, see this branch: https://github.com/tkoolen/drake/tree/tk-mexify

The idea is that you process the mex function inputs based on the argument and return types of the function you want to mexify. The conversion from mxArray to C++ is handled by template <typename T> T fromMex(const mxArray* source), which you specialize for every type you want to be able to handle. The conversion from C++ back to mxArray is handled similarly. I use a variadically templated function mexCallMatlab to recursively bind the first argument of the C++ function to what's returned by fromMex(prhs[0]) until you're left with a zero-argument function (the base case for the recursion), which you can just call.

Here is the relevant code from that branch for your convenience:

#ifndef DRAKE_MEXIFY_H
#define DRAKE_MEXIFY_H

#include <functional>
#include <sstream>
#include <stdexcept>
#include <type_traits>
#include "drakeMexUtil.h"

template <typename T>
T fromMex(const mxArray* mex);

//template <RigidBodyManipulator>
//RigidBody

template<>
RigidBodyManipulator& fromMex(const mxArray* source) {
  return *static_cast<RigidBodyManipulator*>(getDrakeMexPointer(source));
}

template <>
int fromMex(const mxArray* source) {
  return static_cast<int>(mxGetScalar(source));
}

template <>
bool fromMex(const mxArray* source) {
  return mxGetLogicals(source)[0];
}

template <>
std::vector<int>* fromMex(const mxArray* source) {
  return nullptr; // for now
}

template <typename T>
void toMex(const T& source, mxArray* dest[], int nlhs);

template <typename Scalar, int Rows, int Cols>
void toMex(const GradientVar<Scalar, Rows, Cols>& source, mxArray* dest[], int nlhs, bool top_level = true) {
  if (top_level) {
    // check number of output arguments
    if (nlhs > source.maxOrder() + 1) {
      std::ostringstream buf;
      buf << nlhs << " output arguments desired, which is more than the maximum number: " << source.maxOrder() + 1 << ".";
      mexErrMsgTxt(buf.str().c_str());
    }
  }

  if (nlhs != 0)
  {
    // set an output argument
    dest[0] = eigenToMatlab(source.value());

    // recurse
    if (source.hasGradient()) {
      toMex(source.gradient(), &dest[1], nlhs - 1, false);
    }
  }
};

/**
 * mexCallFunction base case
 */
template <typename R>
void mexCallFunction(std::function<R(void)> func, int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
  // call function and convert to mxArray
  toMex(func(), plhs, nlhs);
}

/**
 * mexCallFunction: recursive
 */
template<typename R, typename Arg0, typename ...Args>
void mexCallFunction(std::function<R(Arg0, Args...)> func, int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
  // check number of input arguments
  const int expected_num_args = sizeof...(Args) + 1;
  if (nrhs != expected_num_args) {
    std::ostringstream buf;
    buf << "Expected " << expected_num_args << " arguments, but got " << nrhs << ".";
    mexErrMsgTxt(buf.str().c_str());
  }

  // bind the first argument
  auto partially_applied = [&](Args&& ...tail)
  { return std::move(func)(fromMex<Arg0>(prhs[0]), std::forward<Args>(tail)...); };

  // recursively call mexCallFunction with partially applied function
  mexCallFunction(std::function<R(Args...)>{partially_applied}, nlhs, plhs, nrhs - 1, &prhs[1]);
};


#endif //DRAKE_MEXIFY_H

and the two usages I implemented:

#include "RigidBodyManipulator.h"
#include "mexify.h"

using namespace Eigen;
using namespace std;

void mexFunction(int nlhs, mxArray *plhs[],int nrhs, const mxArray *prhs[]) {
  /*
   * the old implementation was doing a little bit of nontrivial work: inferring the gradient order from nlhs
   * this mexCallFunction implementation hence has a different signature that requires an additional gradient_order input.
   */
  std::function<GradientVar<double, Dynamic, Dynamic>(RigidBodyManipulator&, int)> func = &RigidBodyManipulator::massMatrix<double>;
  mexCallFunction(func, nlhs, plhs, nrhs, prhs);
}

#include "RigidBodyManipulator.h"
#include <mexify.h>

using namespace Eigen;
using namespace std;

void mexFunction( int nlhs, mxArray *plhs[],int nrhs, const mxArray *prhs[] ) {
  /*
   * the old implementation was doing some nontrivial work:
   * 1) creating a vector of ints, passing it into the function, and setting an output argument based on it
   * 2) inferring the gradient order from nlhs
   * this mexCallFunction implementation hence has a different signature and lacks the ability to set the v_or_qdot_indices output
   */
  std::function<GradientVar<double, TWIST_SIZE, Dynamic>(RigidBodyManipulator&, int, int, int, int, bool, std::vector<int>*)> func = &RigidBodyManipulator::geometricJacobian<double>;
  mexCallFunction(func, nlhs, plhs, nrhs, prhs);
}

[RussTedrake]

Very nice. I like it.

At this point, I feel like I understand pretty well how I want the thin wrappers to be implemented, and it would be pretty simple to implement our own swig language to generate the @Class methods given a c++ header file. given Robin’s report, i’m leaning towards the idea that it’s better to do that ourselves than to use the existing swig/matlab bits.

[RussTedrake]

@rdeits -- do you have the .i file for the CoordinateFrame class that I could play with?

[tkoolen]

Glad you like it. I'll try and flesh it out a little more. This does require MSVC 2013 to compile on Windows since variadic templates are unavailable in earlier versions.

[RussTedrake]

i think we’ve agreed to force MSVC >= 2013 at this point.

[rdeits]

@RussTedrake I just started working on a SWIG wrapper file for Drake; I'll pull request the .i file once it's more usable. As a sneak preview, I present the very first working Drake python system:

The .i file for this is pretty short:

%module drake_wrapper

%include <std_except.i>
%include <std_string.i>
%include <windows.i>

%{
#define SWIG_FILE_WITH_INIT
#include <Python.h>
#include "CoordinateFrame.h"
#include "DrakeSystem.h"
#include "Pendulum.h"
%}

%include <typemaps.i>
%include <std_vector.i>
%include <eigen.i>

%template(vectorVectorXd) std::vector<Eigen::VectorXd>;
%template(vectorMatrixXd) std::vector<Eigen::MatrixXd>;
%template(vectorString) std::vector<std::string>;

%eigen_typemaps(Eigen::VectorXd)
%eigen_typemaps(Eigen::MatrixXd)
%eigen_typemaps(Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic>)
%eigen_typemaps(DrakeSystem::VectorXs)

%include "CoordinateFrame.h"
%include "DrakeSystem.h"
%include "Pendulum.h"

[RussTedrake]

woohoo! i love it!

fwiw - i just finished implementing the c++ BotVisualizer class. now i can run the pendulum dynamics in c++, the energy shaping controller in c++, and ddApp
(perhaps the first end-to-end control demo w/ no matlab required)

[RussTedrake]

that .i file is indeed encouragingly simple. exciting.

[rdeits]

By the way, swig recommends having a separate .i file for each module in the python (or whatever language) interface, so we won't have to end up with one mega .i file for all of drake. I'll work on figuring out the right way to structure those files.

[tkoolen]

Wow, that is pretty cool. So eigen.i etc. are already supplied by Swig?

[rdeits]

No, I got the eigen.i file from the internet and then made a few
improvements. I'm keeping the relevant files here:
https://github.com/rdeits/swigmake

On Mon, Aug 24, 2015, 3:38 AM Twan Koolen notifications@github.com wrote:

Wow, that is pretty cool. So eigen.i etc. are already supplied by Swig?

—
Reply to this email directly or view it on GitHub
#1267 (comment)
.

[jaeandersson]

Hi! Joel from CasADi here. As you've noted, there are still some pretty serious memory issues in the module, but we've been making a lot of progress so I'm sure that this and other issues will eventually be resolved. Any help with the module is very much appreciated.
About Octave, SWIG-Octave works very differently from SWIG-MATLAB (despite the fact that SWIG-MATLAB is loosely based on SWIG-Octave). From my understanding, when SWIG-Octave was written, Octave didn't even have class support. This leads to a somewhat unnatural syntax IMO, and although we did get it to work for CasADi, we eventually dropped the idea and focused on SWIG-MATLAB using a fresh module. Since Octave has continued to evolve, I think it's possible to make a new iteration of the SWIG-Octave module. I'm sure the Octave devs would very much welcome this also.

[rdeits]

Hi Joel, it's nice to hear from you! We've still got a few major code decisions to make, and we haven't settled on our favorite way of wrapping our C++ modules. In case you're curious, we've got a few tools and ideas that we've used to try to standardize our mex interfaces somewhat, although we've so far been writing all of the actual mexFunctions by hand.

While swig-matlab uses a single mex file with an internal vector of wrapped objects, we've instead developed DrakeMexPointer which is a wrapper class that contains, at least:

• A raw pointer to memory allocated within the class constructor mexfile
• A deletion function. By convention, deleting a wrapped C++ object means calling the constructor mexfile with the pointer as the sole argument.

We use mexLock and mexUnlock on the constructor mexfile to make sure the memory stays around, and we use mexCallMatlab to create the DrakeMexPointer object from within the mexfile.

This has been pretty stable, although it was a lot of work to get it right. Or, it appeared to be stable, but now we're seeing weird behavior in which the mexLock is perhaps not being respected (??). See: #1274 (comment)

Anyway, we're trying to decide between something similar to our current approach with:

• A mex file for every method
• DrakeMexPointer wrappers for every object, pointing to memory allocated within the mex constructor

vs something more like what swig-matlab does.

Regardless of how we do it, I'd love to end up with something general-purpose and usable beyond our group. Contributing back to swig-matlab seems like a great way to accomplish that, provided we can make it work for our purposes.

[jaeandersson]

OK. There are different ways to generate MATLAB interfaces to C++. The current way is the last iteration in a number of ways that have been attempted to write a MATLAB module to SWIG. And I should say, the most successful approach to date. As you noticed, the current approach uses a single entry point (a single mex function). I think this works very well and since the particular function is chosen by means of a provided index and a switch-block (i.e. constant-time lookup), it is quite efficient as well. An earlier attempt to write a SWIG-MATLAB module was based on "loadlibrary" in MATLAB (which enables multiple entry points), but that approach didn't turn out practical and I think the advantages over mex are negligible. Note that saying that our approach uses a single mex-file is not quite the whole truth: It uses one mex-file per SWIG module and a project can be divided into multiple modules. As of yet, the multiple-module functionality is not however fully implemented.

I don't know your project, but SWIG in general is a great way for improving maintainability of a project, that (like CasADi) has a pretty large API. If for some reason the SWIG-MATLAB module does not fit your needs, note that you can also go half way: Use SWIG to parse your header files and built up a XML-representation of your API. And then implement your own generator from there. We are using this approach successfully to be able to call CasADi from Haskell (which doesn't have a SWIG-module). But I wouldn't go down that road unless you have very good reasons.