I have added a non-linear mixed effects modeling Laplace approximation routine to cppad_py package. One consequence is that one can define an arbitrary Ipopt problem in python and let AD do the derivative calculations. As an example; see https://bradbell.github.io/cppad_py/doc/xsrst/mixed_optimize_fixed_2_py.html