pydrake: Support user-defined dtype for AutoDiffXd and Expression

[EricCousineau-TRI]

EDIT(2020/07/11): Making another top-level summary:

This is one way to achieve reference semantics between C++ and Python for custom dtypes (between Eigen and NumPy), and also some slightly more concrete behavior than what dtype=object offers.
Eric Weister has suggested some ways to override the container (I think), but I had yet to visit that - this will probably be the best way to achieve what I want, I think?
Overall, this has dropped in priority just b/c using dtype=object is sufficiently workable for now, but I would love to revisit this once I have the time.

---

EDIT(2020/06/21): Made upstream issue: pybind/pybind11#2259

---

Per this comment:
pybind/pybind11#1152 (comment)

We can return Eigen::Ref<> in pybind which, if it has a reference return-value-policy and has a parent, will be converted to a numpy.ndarray which is a view to the original data (not a copy of the data). This holds for both immutable and mutable references.

However, custom object data types are presently cast using the effective equivalent of np.array(..., dtype=object), which means the matrix data as a whole will not be aliased.

Resolution is to (a) see if there's a way to do the custom NumPy data types. If that is not feasible, then (b) pydrake should change its semantics to always "mutate" using parent container mutators / return values (which requries extra copying - blech).

Relates #7660, specifically for transmogrification, supporting AutoDiffXd, Symbolic, etc.

---

PR breakdown (as of 2019-01-28), in order:

• numpy: Require NumPy >= 0.15.0 #8608 Incorporating later version of NumPy
• [WIP] Enable NumPy user dtypes #9997 Base support for user-defined dtypes with pybind11
• pydrake autodiff: Use user-defined dtype  #8797 Support for autodiff as user-defined type
• (Some more carving needed)
• [WIP] Use user-defined NumPy dtypes for AutoDiffXd, Expression, etc. #8452 Over-arching feature branch. Contains working prototype code.

Relates:

• Feature requests for Systems Framework Python bindings #7660
• py symbolic: DeprecationWarning not raised by e_xv == e_yv? #10411
• py autodiff: np.sqrt(AutoDiffXd(0.)) encounters divide-by-zero warning #10451
• python symbolic::Formula returns True for matrices (but not scalars)? #8315
• py systems: Ensure usage plays well with caching (invalidation) #7793

[EricCousineau-TRI]

Self notes:

• segfault while casting rational dtype numpy/numpy#5719 - Inadvertently introduces an example of defining a custom NumPy dtype.
• http://www.scipy-lectures.org/advanced/advanced_numpy/ - Docs Mmanu found.

Current plan:

• In a toy pybind11 project, try registering a specific C++ data type in pybind, trying a bit to generalize. Directly including NumPy headers. (Still not sure why they're hidden in pybind11?) Will not directly modify pybind11, if possible.
  • Can overload __eq__ and get different return types.
• See if this can be genearlized / templated.
• Ensure that casting is succinct with what we have with the current functionality (e.g. that the best overloads are respected, and that casting is explicit)
• Try to port to pybind11/eigen.h, or maybe make a separate header, pybind11/numpy_custom.h
  or whatever (seems like the logic is kinda decoupled from Eigen stuff).
• See if it's possible to imply that a non-bool output is first desirable, then see if an explicit option for bool is possible, e.g.

cls.def_ufunc(py::self == py::self);  // CC -> D (non-bool); should be first
cls.def_ufunc(&Class::same_as);  // CC -> B (bool)
# ... In Python
cv == cv  # Returns type of D
np.equal(cv, cv, dtype=bool)  # Returns direct comparison

• See if this change can live next to Add ability to create object matrices pybind/pybind11#1152 to still permit general dtype=object.

Issues:

• (Minor) np.array([...]) will convert to object, even if a custom dtype is registered that's better for it...
  • np.array is defined by multiarraymodule.c, using function _array_fromobject -> PyArray_GetArrayParamsFromObject -> PyArray_DTypeFromObject (recursive inference). There was a point in the code (forgot where) that, if it doesn't find the custom dtype, it just returns NPY_OBJECT :( Could probably be fixed with a wrapping function and a monkey patch, but eww.
  • Actually, numpy.core.test_rational.rational doesn't have this behavior, it works out of the box, Still trying to trace PyArray_DescrConverter and PyArray_DescrCheck

[EricCousineau-TRI]

Posted an update here: #8315 (comment)

Take away is, this seems to be possible, modulo incurring some memory leaks from the NumPy side, per numpy/numpy#10721. There may be some ways around them, but no ways that should affect functionality (as long as we're not depending on the lifetime of our scalars).

[EricCousineau-TRI]

Follow-up response from #8392 (virtualenv stuff):

jwnimmer-tri commented 4 days ago

On the subject of needing numpy fixes (not necessarily virtualenv or python2/3 stuff), I am not sure that I'm following, but my main concern is that in the same way libdrake.so needs to play well with the rest of the system libraries at runtime, because Drake is a library not a program -- pydrake needs to play well with the rest of the python environment at runtime. If we need to use non-default versions of things like numpy, how to we know that other code's use of numpy still works? I think we should go out of our way to be compatible with stable versions of common libraries, and not to require bleeding-edge versions.

jamiesnape commented 4 days ago

I am certainly not in favor of patching packages.

Understood that we do not want bleeding-edge versions, definitely for something as popular for numpy.
The way that pybind11 links to numpy allows for run-time decision to be made, and the only development patches that I am considering are non-critical to functionality (e.g. clean up some memory leaks), and I believe sufficiently small that they are merely implementation changes internal to NumPy, affect no other known behaviors, and a patch simple enough that applying the patch along multiple versions should be tenable.
It is a workflow that I want enough that I'm willing to take responsibility for ensuring this works.

My intent is to do these optional non-API-changing patches until the feature itself reaches upstream numpy, at which point it becomes a bleeding-edge issue that can be resolved via nominal upstream dependency management (which can be alleviated by having an opt-in patch for numpy available).

The specific patch I'd want to implement is to address numpy/numpy#10721.

[jwnimmer-tri]

@EricCousineau-TRI I guess I'm not completely clear on what you're proposing for Drake master as it relates to numpy.

If there are affordances where power users can set their --action_env=PYTHONPATH=... or something to incorporate a more modern numpy at runtime, then I agree that seems like something that you would be able to support on your own, without impacting the rest of the team. (I'd image that would be either no changes to Drake, or at most some #ifdefs in some pybind-related headers.)

If the proposal is to have the default build incorporate such fixes, then I am skeptical that everyone else would be able to ignore what's going on and have the work only fall on your shoulders. If you have a draft commit with what the changes to Drake would look like, that might clear things up for me about what kind of complexity we are talking about.

I don't think any of this prevents you from pushing a bug fix to upstream, but whether and how Drake uses that is the question. Maybe I should wait for an draft of the Drake changes before trying to speculate more?

[jamiesnape]

I agree with @jwnimmer-tri, though my dislike of forks and patches may even be greater. Something would have to be absolutely and undisputedly critical for a patch of something to be justified, other it is opening us up to maintenance and compatibility problems for very moderate gain. I already feel like we are digging ourselves in a hole with pybind11. In most cases, I would say patch upstream and wait or make alternative plans.

[eacousineau]

(Personal acct) Was looking at pybind/pybind11#1731, and stumbled across this:
https://github.com/boostorg/python/blob/b9c0e58/include/boost/python/numpy/ufunc.hpp - will take a deeper peek!

(At a glance, though, it doesn't seem to interface with NumPy UFunc objects; instead, it just writes its own???)

[jwnimmer-tri]

@EricCousineau-TRI a could months ago you changed this from team "kitware" to "manipulation". It seems like a pydrake / bindings question, so should be "kitware"? Or did I misunderstand?

[EricCousineau-TRI]

Er... I don't remember... Happy to switch it back.

[RussTedrake]

Checking in on this one. I'm writing lots of trajectory optimization code for underactuated right now that looks like e.g.:

    for j in range(num_q):  # TODO: use vector form as soon as it's supported
        prog.AddConstraint(qdot_minus[j, i] == qdot_plus[j, i] + h[0, i]*qddot[j])
        prog.AddConstraint(q[j, i+1] == q[j, i] + h[0, i]*qdot_plus[j, i])

I believe that this is the issue that would support adding constraints with the array-type == ? Can you just update me (remind me) of the status here?

[EricCousineau-TRI]

No blockers, just haven't upped it in the priority queue. Want me to schedule this for Q2?

[RussTedrake]

yes, please, if it fits. i think it would add a lot of value.

[jwnimmer-tri]

Note that our minimum required version of numpy is now >= 1.17.4. That might be of help in this issue?