Introducing Magnum Python Bindings

Enter py­bind11

I dis­covered py­bind11 by a lucky ac­ci­dent in early 2018 and im­me­di­ately had to try it. Learn­ing the ba­sics and ex­pos­ing some min­im­al mat­rix/vec­tor math took me about two hours. It was an ex­treme fun and I have to thank all py­bind11 de­velopers for mak­ing it so straight­for­ward to use.

py::class_<Vector3>(m, "Vector3")
    .def_static("x_axis", &Vector3::xAxis, py::arg("length") = 1.0f)
    .def_static("y_axis", &Vector3::yAxis, py::arg("length") = 1.0f)
    .def_static("z_axis", &Vector3::zAxis, py::arg("length") = 1.0f)
    .def(py::init<Float, Float, Float>())
    .def(py::init<Float>())
    .def(py::self == py::self)
    .def(py::self != py::self)
    .def("is_zero", &Vector3::isZero)
    .def("is_normalized", &Vector3::isNormalized)
    .def(py::self += py::self)
    .def(py::self + py::self)
    .def(py::self *= Float{})
    .def(py::self * Float{})
    .def(py::self *= py::self)

How­ever, dif­fer­ent things took a pri­or­ity and so the pro­to­type got shelved un­til it got re­vived again this year. But I learned one main thing — even just the math classes alone were some­thing so use­ful that I kept the built Py­thon mod­ule around and used it from time to time as an en­hanced cal­cu­lat­or. Now, with the mag­num.math mod­ule be­ing al­most com­plete, it’s an every­day tool I use for quick cal­cu­la­tions. Feel free to do the same.

>>> from magnum import *
>>> Matrix3.rotation(Deg(45))
Matrix(0.707107, -0.707107, 0,
    0.707107, 0.707107, 0,
    0, 0, 1)

What Py­thon APIs (and docs) could learn from C++

Every time someone told me they’re us­ing numpy for “do­ing math quickly in Py­thon”, I as­sumed it’s the reas­on­able thing to do — un­til I ac­tu­ally tried to use it. I get that my use case of 4×4 matrices at most might not align well with NumPy’s goals, but the prob­lem is, as far as I know, there’s no full-fea­tured math lib­rary for Py­thon that would give me the whole pack­age1 in­clud­ing Qua­ternions or 2D/3D trans­form­a­tion matrices.

As an ex­cer­cise, for us­ab­il­ity com­par­is­on I tried to ex­press the ro­ta­tion mat­rix shown in the box above in SciPy / NumPy. It took me a good half an hour of star­ing at the docs of scipy.spa­tial.trans­form.Ro­ta­tion un­til I ul­ti­mately de­cided it’s not worth my time. The over­arch­ing prob­lem I have with all those APIs is that it’s not clear at all what types I’m ex­pec­ted to feed to them and provided ex­ample code looks like I’m sup­posed to do half of the cal­cu­la­tions my­self any­way.

>>> from scipy.spatial.transform import Rotation as R
>>> r = R.from_quat([0, 0, np.sin(np.pi/4), np.cos(np.pi/4)])

To avoid the type con­fu­sion, with Mag­num Py­thon bind­ings I de­cided to use strong types where pos­sible — so in­stead of a single dy­nam­ic mat­rix / vec­tor type akin to numpy.ndar­ray, there’s a clear dis­tinc­tion between matrices and vec­tors of dif­fer­ent sizes. So then if you do Matrix4x3() @ Matrix2x4(), docs of Mat­rix4x3.__mat­mul__() will tell you the res­ult is Mat­rix2x3. For NumPy it­self, there’s a pro­pos­al for im­proved type an­nota­tions at numpy/numpy#7370 which would help a lot, but the doc­u­ment­a­tion tools have to make use of that. More on that be­low.

One little thing with big im­pact of the C++ API is strongly-typed angles. You no longer need to re­mem­ber that trig func­tions use ra­di­ans in­tern­ally but HSV col­ors or Open­AL juggles with de­grees in­stead — simply use whatever you please. So Py­thon got the Deg and Rad as well. Py­thon doesn’t have any user-defined lit­er­als (and I’m not aware of any pro­pos­als to add it), how­ever there’s a way to make Py­thon re­cog­nize them. I’m not yet sure if this amount of ma­gic is wise to ap­ply, but I might try it out once.

1.

^ As /u/Ni­hon­Nukite poin­ted out on Red­dit, there’s Pyrr that provides the above miss­ing func­tion­al­ity, fully in­teg­rated with numpy. The only po­ten­tial down­side is that it’s all pure Py­thon, not op­tim­ized nat­ive code.

Hard things are sud­denly easy if you use a dif­fer­ent lan­guage

>>> a = Vector4(1.5, 0.3, -1.0, 1.0)
>>> b = Vector4(7.2, 2.3, 1.1, 0.0)
>>> a.wxy = b.xwz
>>> a
Vector(0, 1.1, -1, 7.2)

If you ever used GLSL or any oth­er shader lan­guage, you prob­ably fell in love with vec­tor swizzles right at the mo­ment you saw them … and then be­came sad after a real­iz­a­tion that such APIs are prac­tic­ally im­pos­sible2 to have in C++. Swizzle op­er­a­tions are nev­er­the­less use­ful and as­sign­ing each com­pon­ent sep­ar­ately would be a pain, so Mag­num provides Math::gath­er() and Math::scat­ter() that al­low you to ex­press the above:

a = Math::scatter<'w', 'x', 'y'>(a, Math::gather<'x', 'w', 'z'>(b));

Verb­ose3 but prac­tic­ally pos­sible. Point is, how­ever, that the above is im­ple­ment­able very eas­ily in Py­thon us­ing __getattr__() and __setattr__() … and a ton of er­ror check­ing on top.

2.

^ GLM does have those, if you en­able GLM_FORCE_SWIZZLE, but do­ing so adds three seconds4 to com­pil­a­tion time of each file that in­cludes GLM head­ers. I’d say that makes swizzles pos­sible in the­ory but such over­head makes them prac­tic­ally use­less.

3.

^ Math func­tions are func­tions and so do not mutate their ar­gu­ments, that’s why the fi­nal self-as­sign­ment. It would of course be bet­ter to be able to write Math::gather<"wxy">(b) or at least Math::gather<'wxy'>(b) but C++ in­sists on the first be­ing im­possible and the second be­ing un­port­able. And cre­at­ing a user-defined lit­er­al just to spe­cify a swizzle seems ex­cess­ive.

4.

^ I did a couple of bench­marks for a yet-to-be-pub­lished art­icle com­par­ing math lib­rary im­ple­ment­a­tions, and this was a shock­er. The only oth­er lib­rary that could come close was Boost.Geo­metry, with two seconds per file.

@overload
def dot(a: Quaternion, b: Quaternion) -> float:
    ...
@overload
def dot(a: Vector2, b: Vector2) -> float:
    ...
def dot(a, b):
    # actual implementation

>>> a = Matrix3.translation((4.0, 2.0))
>>> a
Matrix(1, 0, 4,
       0, 1, 2,
       0, 0, 1)
>>> a.translation = Vector2(5.0, 3.0)
>>> a
Matrix(1, 0, 5,
       0, 1, 3,
       0, 0, 1)

Zero-copy data trans­fer

One very im­port­ant part of Py­thon is the Buf­fer Pro­tocol. It al­lows zero-copy shar­ing of ar­bit­rat­ily shaped data between C and Py­thon — simple tightly-packed lin­ear ar­rays, 2D matrices, or a green chan­nel of a lower right quarter of an im­age flipped up­side down. Hav­ing a full sup­port for the buf­fer pro­tocol was among the reas­ony why Con­tain­ers::StridedAr­rayView went through a ma­jor re­design earli­er this year. This strided ar­ray view is now ex­posed to Py­thon as a con­tain­ers.StridedAr­rayView1D (or Mut­ableStridedAr­rayView1D, and their 2D, 3D and 4D vari­ants) and thanks to the buf­fer pro­tocol it can be seam­lessly con­ver­ted from and to numpy.ar­ray() (and Py­thon’s own memoryview as well). Trans­it­ively that means you can un­leash numpy-based Py­thon al­gorithms dir­ectly on data com­ing out of Im­ageView2D.pixels() and have the modi­fic­a­tions im­me­di­ately re­flec­ted back in C++.

Be­cause, again, hav­ing a spe­cial­ized type with fur­ther re­stric­tions makes the code easi­er to reas­on about, con­tain­ers.Ar­rayView (and its mut­able vari­ant) is ex­posed as well. This one works only with lin­ear tightly packed memory and thus is suit­able for tak­ing views onto bytes or byte­ar­ray, file con­tents and such. Both the strided and lin­ear ar­ray views of course sup­port the full Py­thon sli­cing API. As an ex­ample, here’s how you can read an im­age in Py­thon, pass its con­tents to a Mag­num im­port­er and get the raw pixel data back:

from magnum import trade

def consume_pixels(pixels: np.ndarray):
    ...

importer: trade.AbstractImporter =
    trade.ImporterManager().load_and_instantiate('AnyImageImporter')
with open(filename, 'rb') as f:
    importer.open_data(f.readall())
image: trade.ImageData2D = importer.image2d(0)
# green channel of a lower right quarter of a 256x256 image flipped upside down
consume_pixels(image.pixels[128:128,::-1,1:2])

Just one ques­tion left — who owns the memory here, then? To an­swer that, let’s dive in­to Py­thon’s ref­er­ence count­ing.

Ref­er­en­ce count­ing

In C++, views are one of the more dan­ger­ous con­tain­ers, as they ref­er­ence data owned by some­thing else. There you’re ex­pec­ted to en­sure the data own­er is kept in scope for at least as long as the view on it. A sim­il­ar thing is with oth­er types — for ex­ample, a GL::Mesh may ref­er­ence a bunch of GL::Buf­fers, or a Trade::Ab­strac­tIm­port­er loaded from a plu­gin needs its plu­gin man­ager to be alive to keep the plu­gin lib­rary loaded.

How­ever, im­pos­ing sim­il­ar con­straints on Py­thon users would be dar­ing too much, so all ex­posed Mag­num types that refer to ex­tern­al data im­ple­ment ref­er­ence count­ing un­der the hood. The des­ig­nated way of do­ing this with py­bind11 is wrap­ping all your everything with std::shared_ptr. On the oth­er hand, Mag­num is free of any shared point­ers by design, and adding them back just to make Py­thon happy would make every­one else angry in ex­change. What Mag­num does in­stead is ex­tend­ing the so-called hold­er type in py­bind11 (which doesn’t have to be std::shared_ptr; std::unique_ptr or a cus­tom point­er types is fine as well) and stor­ing ref­er­ences to in­stance de­pend­en­cies in­side it.

The straight­for­ward way of do­ing this would be to take GL::Mesh, sub­class it in­to a PyMesh, store buf­fer ref­er­ences in­side it and then ex­pose PyMesh as gl.Mesh in­stead. But com­pared to the hold­er type ap­proach this has a ser­i­ous dis­ad­vant­age where every API that works with meshes would sud­denly need to work with PyMesh in­stead and that’s not al­ways pos­sible.

For test­ing and de­bug­ging pur­poses, ref­er­ences to memory own­ers or oth­er data are al­ways ex­posed through the API — see for ex­ample Im­ageView2D.own­er or gl.Mesh.buf­fers.

view = containers.ArrayView(data)
while view:
    # Comment, ignore until EOL
    if view[0] == '#': while view and view[0] != '\n': view = view[1:]
    # Vertex / face
    elif view[0] == 'v': view = self.parse_vertex(view)
    elif view[0] == 'f': view = self.parse_face(view)
    ...

The less-than-great as­pects of py­bind11

.def("__getitem__", [](const T& self, std::size_t i) -> typename T::Type {
    if(i >= T::Size) throw pybind11::index_error{};
    return self[i];
})

struct buffer_info {
    void *ptr;
    ssize_t itemsize;
    std::string format;
    ssize_t ndim;
    std::vector<ssize_t> shape;
    std::vector<ssize_t> strides;
};

Every­one “just uses Sphinx”. You?

The ob­vi­ous first choice when it comes to doc­u­ment­ing Py­thon code is to use Sphinx — everything in­clud­ing the stand­ard lib­rary uses it and I don’t even re­mem­ber see­ing a single Py­thon lib­rary that doesn’t. How­ever, if you clicked on any of the above doc links, you prob­ably real­ized that … no, Mag­num is not us­ing it.

Ever since the doc­u­ment­a­tion search got in­tro­duced early last year, many de­velopers quickly be­came ad­dicted used to it. Whip­ping up some Sphinx docs, where both search per­form­ance and res­ult rel­ev­ance is ex­tremely un­der­whelm­ing, would be ef­fect­ively un­do­ing all us­ab­il­ity pro­gress Mag­num made un­til now, so the only op­tion was to bring the search to Py­thon as well.

While at it, I made the doc gen­er­at­or aware of all kinds of type an­nota­tions, prop­erly cross­link­ing everything to cor­res­pond­ing type defin­i­tions. And not just loc­al types — sim­il­arly to Doxy­gen’s tag­files, Sphinx has In­ter­sphinx, so link­ing to 3rd party lib­rary docs (such as NumPy) or even the stand­ard lib­rary is pos­sible as well. Con­versely, the m.css Py­thon doc gen­er­at­or ex­ports an In­ter­sphinx in­vent­ory file, so no mat­ter wheth­er you use vanilla Sphinx or m.css, you can link to m.css-gen­er­ated Py­thon docs from your own doc­u­ment­a­tion as well.

If you want to try it out on your pro­ject, head over to the m.css web­site for a brief in­tro­duc­tion. Com­pared to Sphinx or Doxy­gen it be­haves more like Doxy­gen, as it im­pli­citly browses your mod­ule hier­archy and gen­er­ates a ded­ic­ated page for each class and mod­ule. The way Sphinx does it was in­ter­est­ing for me at first, but over the time I real­ized it needs quite a lot of ef­fort from de­veloper side to or­gan­ize well — and from the doc­u­ment­a­tion read­er side, it can lead to things be­ing harder to find than they should be (for ex­ample, docs for str.splitlines() are bur­ied some­where in the middle of a kilo­met­er-long page doc­u­ment­ing all built­in types).

The doc gen­er­at­or re­sembles Sphinx, but I de­cided to ex­per­i­ment with a clean slate first in­stead of mak­ing it 100% com­pat­ible — some design de­cisions in Sphinx it­self are his­tor­ic­al (such as type an­nota­tions in the doc block in­stead of in the code it­self) and it didn’t make sense for me to port those over. Well, at least for now, a full Sphinx com­pat­ib­il­ity is not com­pletely out of ques­tion.

What’s next?

The work is far from be­ing done — apart from ex­pos­ing APIs that are not ex­posed yet (which is just a routine work, mostly), I’m still not quite sat­is­fied with the per­form­a­ce of bound types, so on my roadmap is try­ing to ex­pose a ba­sic type us­ing pure C Py­thon APIs the most ef­fi­cient way pos­sible and then com­par­ing how long does it take to in­stan­ti­ate that type and call meth­ods on it. One of the things to try is a vec­tor­call call pro­tocol that’s new in Py­thon 3.8 (PEP590) and the re­search couldn’t be com­plete without also try­ing a sim­il­ar thing in Mi­croPy­thon.

~ ~ ~