Adds a NumPy backend so opty can be used without Cython (and thus C compilation).

[moorepants]

This is a first shot at fixing #367 but I quickly hit a limitation of trying to have the loop evaluation purely as a NumPy vectorized operation. See the following error:

(opty-dev) moorepants@nandi:opty(master)$ py.test opty/tests/test_direct_collocation.py 
==================================================================================== test session starts =====================================================================================
platform linux -- Python 3.13.1, pytest-8.3.4, pluggy-1.5.0
rootdir: /home/moorepants/src/opty
plugins: cov-6.0.0
collected 64 items                                                                                                                                                                           

opty/tests/test_direct_collocation.py ...........F....................................................                                                                                 [100%]

========================================================================================== FAILURES ==========================================================================================
__________________________________________________________ TestConstraintCollocator.test_gen_multi_arg_con_jac_func_backward_euler ___________________________________________________________

self = <opty.tests.test_direct_collocation.TestConstraintCollocator object at 0x70de03402a80>

    def test_gen_multi_arg_con_jac_func_backward_euler(self):
    
        self.collocator._gen_multi_arg_con_jac_func()
        self.numpy_collocator._gen_multi_arg_con_jac_func()
    
        # Make sure the parameters are in the correct order.
        constant_values = \
            np.array([self.constant_values[self.constant_symbols.index(c)]
                      for c in self.collocator.parameters])
    
        # TODO : Once there are more than one specified, they will need to
        # be put in the correct order too.
    
        jac_vals = self.collocator._multi_arg_con_jac_func(self.state_values,
                                                           self.specified_values,
                                                           constant_values,
                                                           self.interval_value)
>       numpy_jac_vals = self.numpy_collocator._multi_arg_con_jac_func(
            self.state_values, self.specified_values, constant_values,
            self.interval_value)

opty/tests/test_direct_collocation.py:332: 
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 
opty/direct_collocation.py:2041: in constraints_jacobian
    non_zero_derivatives = eval_partials(*args)
opty/utils.py:594: in transpose
    return np.moveaxis(eval_single_mat(*array_num_args), -1, 0)
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 

xi = array([2., 3., 4.]), vi = array([6., 7., 8.]), xp = array([1., 2., 3.]), vp = array([5., 6., 7.]), fi = array([2., 2., 2.]), m = array([1., 1., 1.]), c = array([2., 2., 2.])
k = array([3., 3., 3.]), h_opty = array([0.01, 0.01, 0.01])

    def _lambdifygenerated(xi, vi, xp, vp, fi, m, c, k, h_opty):
>       return array([[h_opty**(-1.0), -1, -1/h_opty, 0, 0], [k, c + m/h_opty, 0, -m/h_opty, xi]])
E       ValueError: setting an array element with a sequence. The requested array has an inhomogeneous shape after 2 dimensions. The detected shape was (2, 5) + inhomogeneous part.

<lambdifygenerated-9>:2: ValueError
================================================================================== short test summary info ===================================================================================
FAILED opty/tests/test_direct_collocation.py::TestConstraintCollocator::test_gen_multi_arg_con_jac_func_backward_euler - ValueError: setting an array element with a sequence. The requested array has an inhomogeneous shape after 2 dimensions. The detected shape was (2, 5) + inhomogeneous part.
========================================================================== 1 failed, 63 passed in 151.01s (0:02:31) ==========================================================================

and then the SymPy issues that explain why this can't be done: sympy/sympy#27632

I will change the code to function in the next commits, but the performance will likely suffer.

[moorepants]

Timing comparison for the one legged time trial example:

NumPy backend:

In [2]: %timeit problem.jacobian(initial_guess)
90.3 ms ± 2.07 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

In [3]: %timeit problem.constraints(initial_guess)
38.2 ms ± 1.97 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

Cython backend:

In [2]: %timeit problem.jacobian(initial_guess)
748 μs ± 26.2 μs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)

In [3]: %timeit problem.constraints(initial_guess)
306 μs ± 3.09 μs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)

[Peter230655]

Numpy looks much faster or am I reading it wrong?

[moorepants]

NumPy is 120x slower for this example.

[Peter230655]

I mixed up millisec and microsec!☹️

[moorepants]

Do you have any thoughts on this new api? backend='cython'|'numpy' ?

[moorepants]

I still get a warning for a pow() when evaluating. I'm not sure what number is causing that.

[moorepants]

opty/tests/test_utils.py::test_ufuncify_matrix
  /home/moorepants/src/opty/opty/tests/test_utils.py:283: RuntimeWarning: invalid value encountered in power
    result[:, 0, 0] = a_vals**2 * np.cos(np.pi*b_vals)**c_vals

opty/tests/test_utils.py::test_ufuncify_matrix
  <lambdifygenerated-25>:4: RuntimeWarning: invalid value encountered in scalar power
    return array([[x0*cos(pi*b)**_Dummy_250, _Dummy_250**4 + tan(b)/sin(x1)], [-sqrt(_Dummy_250) + b**2 + x0, x1*(_Dummy_250 + x1)*sin(b)/a]])

[Peter230655]

Do you have any thoughts on this new api? backend='cython'|'numpy' ?

I assume for a computer idiot like me, numpy is better: one less 'thing' to mess around with.
(When I first installed this stuff, Timo wrote detailed instructions, which my son could follow and set it up)
As I understand, numpy uses C and FORTRAN so it should be fast - but I as understood your disucssion with Aaron, you cannot use the full speed of numpy?

[moorepants]

I am implementing this so you and others can use opty without having to install a C compiler and ensure the compiler is working properly with Python (an issue that you and others have faced). So, yes I'm trying to make the software simpler to use.

Yes, the conversation with Aaron points to the fact that I can't harness the full speed of NumPy when using SymPy's lambdify for this specific problem. NumPy will actually run these computations faster than what is shown here, but currently not via SymPy's lambdify code generation. The fact that I have to have a loop through all nodes in Python code, is the reason it is slow.

But, I'm not really asking you about either of those points. I'm curious if you think how I've added this feature is appropriate for opty. Did I design the API well? Should I use different variable names? Should I have different documentation? i.e can your review these code changes because it is good that we feel pretty confident that we wouldn't need to change any of the API later (causing backwards incompatibilities).

[moorepants]

Also, do I have enough tests? Does it cover all cases? etc.

[Peter230655]

I am implementing this so you and others can use opty without having to install a C compiler and ensure the compiler is working properly with Python (an issue that you and others have faced). So, yes I'm trying to make the software simpler to use.

Yes, the conversation with Aaron points to the fact that I can't harness the full speed of NumPy when using SymPy's lambdify for this specific problem. NumPy will actually run these computations faster than what is shown here, but currently not via SymPy's lambdify code generation. The fact that I have to have a loop through all nodes in Python code, is the reason it is slow.

But, I'm not really asking you about either of those points. I'm curious if you think how I've added this feature is appropriate for opty. Did I design the API well? Should I use different variable names? Should I have different documentation? i.e can your review these code changes because it is good that we feel pretty confident that we wouldn't need to change any of the API later (causing backwards incompatibilities).

So that I understand correctly:
I will look at the docstrings of the changes you made and comment? (You did not put anything in README.rst as fa as I could tell.)

[moorepants]

I'm asking you to look at all the changes I made in this PR and give specific feedback on the changes.

[Peter230655]

I'm asking you to look at all the changes I made in this PR and give specific feedback on the changes.

Clear! Already started.

In line 586 you write function
In line 587 you write Function. Maybe better to also write function?

in direct_colocation.py
line 785: backend
Line 786: Backend. Maybe better backend here, too - after all this is its name.
(I like the term backend!)

Dumb question:
Would backend (and cse) not also be optional kwargs of Problem?
The user 'communicates' with opty through Problem.

[moorepants]

For opty, I do not think there is a reason a user would ever not want cse, so it it has never been an option for them to say anything about that. Before the addition of this function, there was no choice at all to not have cse. I actually only added it to lambdify_matrix() because I was hitting some bugs with lambdify and needed to disable it and enable it to debug. A user should never run lambdify_matrix or ufuncify_matrix and I purposely don't expose them in the documentation. So, I guess I should just remove cse kwarg from lambdify_matrix.

Note that you can press the blue "+" symbol on any line in the files changes tab to add a comment directly at that line.

[moorepants]

Also, thanks for the very helpful feedback. Exactly what I was hoping for!

[moorepants]

The capitalization here is because "backend" is the first word in the sentence. Just above it is lowercase because that is the exact case-sensitive kwarg. @Peter230655 can you clarify you comment on this?

[Peter230655]

The capitalization here is because "backend" is the first word in the sentence. Just above it is lowercase because that is the exact case-sensitive kwarg. @Peter230655 can you clarify you comment on this?

What I meant was this: backend, with small b is the name of the kwarg. Hence, no matter where it stands in a sentence its name should be spelled the same way. But a VERY minor point indeed!

[moorepants]

I'm using "Backend" in the sentence as normal English word, like "The backend of the government's computer system is very complex."

[Peter230655]

For opty, I do not think there is a reason a user would ever not want cse, so it it has never been an option for them to say anything about that. Before the addition of this function, there was no choice at all to not have cse. I actually only added it to lambdify_matrix() because I was hitting some bugs with lambdify and needed to disable it and enable it to debug. A user should never run lambdify_matrix or ufuncify_matrix and I purposely don't expose them in the documentation. So, I guess I should just remove cse kwarg from lambdify_matrix.

Note that you can press the blue "+" symbol on any line in the files changes tab to add a comment directly at that line.

I totally agree regardind cse=True always: You may remember that I wotze some sympy mechanics simulations with rather large eoms. Only after I learned about cse=True to be used in lambdify would solve_ivp get a result in good time.

I did not know about this blue "+". Let me try it.

[moorepants]

I think I'm good with this, so merging. We can tweak before a next release if needed.