Switched to Cython memoryviews

examples-gallery/plot_parallel_park.py

@@ -8,7 +8,7 @@
 **Constants**
 
 - m : car mass, [kg]
-- I : car yaw moment of inertia, [kg m^2]
+- Izz : car yaw moment of inertia, [kg m^2]
 - a : distance from front axle to mass center, [m]
 - b : distance from rear axle to mass center, [m]
 
@@ -40,7 +40,7 @@
 
 # %%
 # Generate the nonholonomic equations of motion of the system.
-m, I, a, b = sm.symbols('m, I, a, b', real=True)
+m, Izz, a, b = sm.symbols('m, Izz, a, b', real=True)
 x, y, vx, vy = me.dynamicsymbols('x, y, v_x, v_y', real=True)
 theta, omega = me.dynamicsymbols('theta, omega', real=True)
 delta, beta = me.dynamicsymbols('delta, beta', real=True)
@@ -76,9 +76,9 @@
     Pf.vel(N).dot(B.y),
 ]
 
-IA = me.inertia(A, 0, 0, I)
+IA = me.inertia(A, 0, 0, Izz)
 car = me.RigidBody('A', Ao, A, m, (IA, Ao))
-IB = me.inertia(B, 0, 0, I/32)
+IB = me.inertia(B, 0, 0, Izz/32)
 wheel = me.RigidBody('B', Pf, B, m/6, (IB, Pf))
 
 propulsion = (Pr, F*A.x)
@@ -112,7 +112,7 @@
 # %%
 # Provide some reasonably realistic values for the constants.
 par_map = {
-    I: 1/12*1200*(2**2 + 3**2),
+    Izz: 1/12*1200*(2**2 + 3**2),
     m: 1200,
     a: 1.5,
     b: 1.5,

examples-gallery/plot_pendulum_swing_up_fixed_duration.py

@@ -23,21 +23,21 @@
 
 # %%
 # Specify the symbolic equations of motion.
-I, m, g, d, t = sm.symbols('I, m, g, d, t')
+Izz, m, g, d, t = sm.symbols('Izz, m, g, d, t')
 theta, omega, T = sm.symbols('theta, omega, T', cls=sm.Function)
 
 state_symbols = (theta(t), omega(t))
-constant_symbols = (I, m, g, d)
+constant_symbols = (Izz, m, g, d)
 specified_symbols = (T(t),)
 
 eom = sm.Matrix([theta(t).diff() - omega(t),
-                 I*omega(t).diff() + m*g*d*sm.sin(theta(t)) - T(t)])
+                 Izz*omega(t).diff() + m*g*d*sm.sin(theta(t)) - T(t)])
 sm.pprint(eom)
 
 # %%
 # Specify the known system parameters.
 par_map = {
-    I: 1.0,
+    Izz: 1.0,
     m: 1.0,
     g: 9.81,
     d: 1.0,
@@ -51,7 +51,8 @@
 obj, obj_grad = create_objective_function(obj_func, state_symbols,
                                           specified_symbols, tuple(),
                                           num_nodes,
-                                          node_time_interval=interval_value)
+                                          node_time_interval=interval_value,
+                                          time_symbol=t)
 
 # %%
 # Specify the symbolic instance constraints, i.e. initial and end conditions,
@@ -71,11 +72,11 @@
 
 # %%
 # Create an optimization problem.
-prob = Problem(obj, obj_grad, eom, state_symbols,
-               num_nodes, interval_value,
+prob = Problem(obj, obj_grad, eom, state_symbols, num_nodes, interval_value,
                known_parameter_map=par_map,
                instance_constraints=instance_constraints,
-               bounds=bounds)
+               bounds=bounds,
+               time_symbol=t)
 
 # %%
 # Use a random positive initial guess.

opty/tests/test_direct_collocation.py

@@ -19,7 +19,7 @@ def test_pendulum():
     interval_value = duration / (num_nodes - 1)
 
     # Symbolic equations of motion
-    I, m, g, d, t = sym.symbols('I, m, g, d, t')
+    I, m, g, d, t = sym.symbols('Ix, m, g, d, t')
     theta, omega, T = sym.symbols('theta, omega, T', cls=sym.Function)
 
     state_symbols = (theta(t), omega(t))
@@ -818,7 +818,7 @@ class TestConstraintCollocatorInstanceConstraints():
 
     def setup_method(self):
 
-        I, m, g, d, t = sym.symbols('I, m, g, d, t')
+        I, m, g, d, t = sym.symbols('Ix, m, g, d, t')
         theta, omega, T = [f(t) for f in sym.symbols('theta, omega, T',
                                                      cls=sym.Function)]
 

opty/tests/test_utils.py

@@ -65,7 +65,7 @@ def setup_method(self):
         self.q = len(self.input_symbols)
         self.r = len(self.unknown_symbols)
         self.N = 20
-        
+
         self.x_vals = np.random.random(self.N)
         self.v_vals = np.random.random(self.N)
         self.f1_vals = np.random.random(self.N)
@@ -109,7 +109,7 @@ def test_backward_all(self):
         obj_expr = (
             sym.Integral(self.x ** 2 + self.m ** 2, self.t) +
             sym.Integral(self.c ** 2 * self.f2 ** 2 , self.t) +
-            sym.sin(self.k) ** 2            
+            sym.sin(self.k) ** 2
         )
         obj, obj_grad = utils.create_objective_function(
             obj_expr, self.state_symbols, self.input_symbols,
@@ -163,7 +163,7 @@ def test_midpoint_all(self):
         obj_expr = (
             sym.Integral(self.x ** 2 + self.m ** 2, self.t) +
             sym.Integral(self.c ** 2 * self.f2 ** 2 , self.t) +
-            sym.sin(self.k) ** 2            
+            sym.sin(self.k) ** 2
         )
         x_mid = (self.x_vals[1:] + self.x_vals[:-1]) / 2
         f2_mid = (self.f2_vals[1:] + self.f2_vals[:-1]) / 2

opty/utils.py

@@ -252,9 +252,10 @@ def parse_free(free, n, q, N, variable_duration=False):
         return free_states, free_specified, free_constants
 
 
-def create_objective_function(
-    objective, state_symbols, input_symbols, unknown_symbols,
-    N, node_time_interval=1.0, integration_method="backward euler"):
+def create_objective_function(objective, state_symbols, input_symbols,
+                              unknown_symbols, N, node_time_interval=1.0,
+                              integration_method="backward euler",
+                              time_symbol=me.dynamicsymbols._t):
     """Creates function to evaluate the objective and objective gradient.
 
     Parameters
@@ -278,6 +279,8 @@ def create_objective_function(
     integration_method : str, optional
         The method used to integrate the system. The default is "backward
         euler".
+    time_symbol : sympy.Symbol
+        If not supplied, ``sympy.physics.mechanics.dynamicsymbols._t`` is used.
 
     """
     def lambdify_function(expr, multiplication_array, take_sum):
@@ -298,7 +301,7 @@ def parse_expr(expr, in_integral=False):
         if isinstance(expr, sm.Integral):
             if in_integral:
                 raise NotImplementedError("Nested integrals are not supported.")
-            if expr.limits != ((me.dynamicsymbols._t,),):
+            if expr.limits != ((time_symbol,),):
                 raise NotImplementedError(
                     "Only indefinite integrals of time are supported.")
             return int_placeholder(parse_expr(expr.function, True))
@@ -413,18 +416,24 @@ def sort_sympy(seq):
 cimport cython
 
 cdef extern from "{file_prefix}_h.h"{head_gil}:
-    void {routine_name}(double matrix[{matrix_output_size}], {input_args})
+    void {routine_name}(double matrix[{matrix_output_size}],
+{input_args})
 
 @cython.boundscheck(False)
 @cython.wraparound(False)
-def {routine_name}_loop(np.ndarray[np.double_t, ndim=2] matrix, {numpy_typed_input_args}):
+def {routine_name}_loop(matrix,
+{numpy_typed_input_args}):
+
+    cdef double[:, ::1] matrix_memview = matrix
+    {memory_views}
 
     cdef int n = matrix.shape[0]
 
     cdef int i
 
     for i in {loop_sig}:
-        {routine_name}(&matrix[i, 0], {indexed_input_args})
+        {routine_name}(&matrix_memview[i, 0],
+{indexed_input_args})
 
     return matrix.reshape(n, {num_rows}, {num_cols})
 """
@@ -604,31 +613,41 @@ def ufuncify_matrix(args, expr, const=None, tmp_dir=None, parallel=False,
     d['eval_code'] = '    ' + '\n    '.join((sub_expr_code + '\n' +
                                              matrix_code).split('\n'))
 
-    c_indent = len('void {routine_name}('.format(**d))
+    c_indent = len('    void {routine_name}('.format(**d))
     c_arg_spacer = ',\n' + ' ' * c_indent
 
     input_args = ['double {}'.format(sm.ccode(a)) for a in args]
-    d['input_args'] = c_arg_spacer.join(input_args)
+    d['input_args'] = ' '*c_indent + c_arg_spacer.join(input_args)
 
     cython_input_args = []
     indexed_input_args = []
+    memory_views = []
     for a in args:
         if const is not None and a in const:
             typ = 'double'
             idexy = '{}'
+            cython_input_args.append('{} {}'.format(typ, sm.ccode(a)))
         else:
-            typ = 'np.ndarray[np.double_t, ndim=1]'
-            idexy = '{}[i]'
+            idexy = '{}_memview[i]'
+            memview = 'cdef double[::1] {}_memview = {}'
+            memory_views.append(memview.format(sm.ccode(a), sm.ccode(a)))
+            cython_input_args.append('{}'.format(sm.ccode(a)))
 
-        cython_input_args.append('{} {}'.format(typ, sm.ccode(a)))
         indexed_input_args.append(idexy.format(sm.ccode(a)))
 
     cython_indent = len('def {routine_name}_loop('.format(**d))
-    cython_arg_spacer = ',\n' + ' ' * cython_indent
+    cython_arg_spacer = ',\n' + ' '*cython_indent
+
+    loop_indent = len('        {routine_name}('.format(**d))
+    loop_spacer = ',\n' + ' '*loop_indent
+
+    d['numpy_typed_input_args'] = (' '*cython_indent +
+                                   cython_arg_spacer.join(cython_input_args))
 
-    d['numpy_typed_input_args'] = cython_arg_spacer.join(cython_input_args)
+    d['indexed_input_args'] = (' '*loop_indent +
+                               loop_spacer.join(indexed_input_args))
 
-    d['indexed_input_args'] = ',\n'.join(indexed_input_args)
+    d['memory_views'] = '\n    '.join(memory_views)
 
     files = {}
     files[d['file_prefix'] + '_c.c'] = _c_template.format(**d)