Problems with sympy.pi in the equations of motion

[Peter230655]

This morning I ran into a problem using opty for a simulation (Example 10.59 from Betts' book)
With the help from @tjstienstra I found the problem was that I used sympy.pi in the equations of motion. Replacing it with numpy.pi solved the problem.
Here I think, one can see that opty objects to sympy.pi:

in the simulation below, I added a factor sympy.cos(symypy.pi) / sympy.cos(numpy.pi) and sympy.pi / numpy.py, in lines 37, 38.
For some reason, it accepts sympy.cos(sympy.pi), but not sympy.pi alone.

NB:

1. As I recall, I had this problem before, but I do not think, I raised an issue then.
2. Of course no big deal if one knows it, but maybe this could be mentioned in the docs?

"""
Parameter Identification: Betts & Huffman 2003
==============================================

This is the problem presented in section 7 of:

Betts, John T. and Huffman, William P. "Large Scale Parameter Estimation
Using Sparse Nonlinear Programming Methods". SIAM J. Optim., Vol 14, No. 1,
pp. 223-244, 2003.

"""

from collections import OrderedDict

import numpy as np
import sympy as sym
import matplotlib.pyplot as plt
from opty.direct_collocation import Problem

duration = 1.0
num_nodes = 100
interval = duration / (num_nodes - 1)

# %%
# Symbolic equations of motion
mu, p, t = sym.symbols('mu, p, t')
y1, y2, T = sym.symbols('y1, y2, T', cls=sym.Function)

state_symbols = (y1(t), y2(t))
constant_symbols = (mu, p)

# %%
# I had to use a little "trick" here because time was explicit in the eoms,
# i.e. setting T(t) to a function of time and then pass in time as a known
# trajectory in the problem.
eom = sym.Matrix([y1(t).diff(t) - y2(t)*sym.cos(sym.pi)/sym.cos(np.pi),
                  y2(t).diff(t) - mu**2 * y1(t) + (mu**2 + p**2) *
                  sym.sin(p * T(t))*sym.pi/np.pi])

# %%
# Specify the known system parameters.
par_map = OrderedDict()
par_map[mu] = 60.0

# %%
# Generate data
time = np.linspace(0.0, 1.0, num_nodes)
y1_m = np.sin(np.pi * time) + np.random.normal(scale=0.05, size=len(time))
y2_m = np.pi * np.cos(np.pi * time) + np.random.normal(scale=0.05,
                                                       size=len(time))


# %%
# Specify the objective function and it's gradient. I'm only fitting to y1,
# but they may have fit to both states in the paper.
def obj(free):
    return interval * np.sum((y1_m - free[:num_nodes])**2)


def obj_grad(free):
    grad = np.zeros_like(free)
    grad[:num_nodes] = 2.0 * interval * (free[:num_nodes] - y1_m)
    return grad


# %%
# Specify the symbolic instance constraints, i.e. initial and end
# conditions.
instance_constraints = (y1(0.0), y2(0.0) - np.pi)

# %%
# Create an optimization problem.
prob = Problem(obj, obj_grad,
               eom, state_symbols,
               num_nodes, interval,
               known_parameter_map=par_map,
               known_trajectory_map={T(t): time},
               instance_constraints=instance_constraints,
               time_symbol=t,
               integration_method='midpoint')

# %%
# Use a random positive initial guess.
initial_guess = np.random.randn(prob.num_free)

# %%
# Find the optimal solution.
solution, info = prob.solve(initial_guess)
print(info['status_msg'])

# %%
# Print results.
known_msg = "Known value of p = {}".format(np.pi)
identified_msg = "Identified value of p = {}".format(solution[-1])
divider = '=' * max(len(known_msg), len(identified_msg))

print(divider)
print(known_msg)
print(identified_msg)
print(divider)

# %%
# Plot results
fig_y1, axes_y1 = plt.subplots(3, 1)

legend = ['measured', 'initial guess', 'direct collocation solution']

axes_y1[0].plot(time, y1_m, '.k',
                time, initial_guess[:num_nodes], '.b',
                time, solution[:num_nodes], '.r')
axes_y1[0].set_xlabel('Time [s]')
axes_y1[0].set_ylabel('y1')
axes_y1[0].legend(legend)

axes_y1[1].set_title('Initial Guess Constraint Violations')
axes_y1[1].plot(prob.con(initial_guess)[:num_nodes - 1])
axes_y1[2].set_title('Solution Constraint Violations')
axes_y1[2].plot(prob.con(solution)[:num_nodes - 1])

plt.tight_layout()

fig_y2, axes_y2 = plt.subplots(3, 1)

axes_y2[0].plot(time, y2_m, '.k',
                time, initial_guess[num_nodes:-1], '.b',
                time, solution[num_nodes:-1], '.r')
axes_y2[0].set_xlabel('Time [s]')
axes_y2[0].set_ylabel('y2')
axes_y2[0].legend(legend)

axes_y2[1].set_title('Initial Guess Constraint Violations')
axes_y2[1].plot(prob.con(initial_guess)[num_nodes - 1:])
axes_y2[2].set_title('Solution Constraint Violations')
axes_y2[2].plot(prob.con(solution)[num_nodes - 1:])

plt.tight_layout()

plt.show()

[tjstienstra]

Here is a minimal reproducer of the issue:

import sympy as sm
import numpy as np
from opty.utils import ufuncify_matrix

x, y = sm.symbols("x y")
result, x_vals, y_vals = np.empty((1, 2)), np.array([1.0]), np.array([2.0])
expected = np.array([[1.0, 6.283185307179586]])
f_np = ufuncify_matrix((x, y), sm.Matrix([x, np.pi * y]))  # works fine
f_np(result, x_vals, y_vals)
np.testing.assert_allclose(result, expected)
f_sm = ufuncify_matrix((x, y), sm.Matrix([x, sm.pi * y]))  # crashes due to a compilation error
f_sm(x_vals, y_vals, result)
np.testing.assert_allclose(result, expected)

[moorepants]

When I run Timo's example on Linux I get this error:

Traceback (most recent call last):
  File "/home/moorepants/src/opty/tmp.py", line 12, in <module>
    f_sm(x_vals, y_vals, result)
    ~~~~^^^^^^^^^^^^^^^^^^^^^^^^
  File "ufuncify_matrix_1.pyx", line 17, in ufuncify_matrix_1.eval_matrix_loop
ValueError: Buffer has wrong number of dimensions (expected 2, got 1)

[Peter230655]

I vaguely recall, that my example above (or rather something similar) did not give an error witth Linux, only with windows - but maybe my memory is wrong.

[Peter230655]

When I run Timo's example, I get this error:

---------------------------------------------------------------------------
ModuleNotFoundError                       Traceback (most recent call last)
File c:\users\peter\github_repos\opty\opty\utils.py:744, in ufuncify_matrix(args, expr, const, tmp_dir, parallel, show_compile_output)
    [743](file:///C:/users/peter/github_repos/opty/opty/utils.py:743) try:
--> [744](file:///C:/users/peter/github_repos/opty/opty/utils.py:744)     cython_module = importlib.import_module(d['file_prefix'])
    [745](file:///C:/users/peter/github_repos/opty/opty/utils.py:745) except ImportError as error:

File c:\Users\Peter\anaconda3\envs\opty-dev\Lib\importlib\__init__.py:90, in import_module(name, package)
     [89](file:///C:/Users/Peter/anaconda3/envs/opty-dev/Lib/importlib/__init__.py:89)         level += 1
---> [90](file:///C:/Users/Peter/anaconda3/envs/opty-dev/Lib/importlib/__init__.py:90) return _bootstrap._gcd_import(name[level:], package, level)

File <frozen importlib._bootstrap>:1387, in _gcd_import(name, package, level)

File <frozen importlib._bootstrap>:1360, in _find_and_load(name, import_)

File <frozen importlib._bootstrap>:1324, in _find_and_load_unlocked(name, import_)

ModuleNotFoundError: No module named 'ufuncify_matrix_1'

The above exception was the direct cause of the following exception:

ImportError                               Traceback (most recent call last)
File Untitled-1:12
     [10](untitled-1:10) f_np(result, x_vals, y_vals)
     [11](untitled-1:11) np.testing.assert_allclose(result, expected)
---> [12](untitled-1:12) f_sm = ufuncify_matrix((x, y), sm.Matrix([x, sm.pi * y]))  # crashes due to a compilation error
     [13](untitled-1:13) f_sm(x_vals, y_vals, result)
     [14](untitled-1:14) np.testing.assert_allclose(result, expected)

File c:\users\peter\github_repos\opty\opty\utils.py:749, in ufuncify_matrix(args, expr, const, tmp_dir, parallel, show_compile_output)
    [745](file:///C:/users/peter/github_repos/opty/opty/utils.py:745)     except ImportError as error:
    [746](file:///C:/users/peter/github_repos/opty/opty/utils.py:746)         msg = ('Unable to import the compiled Cython module {}, '
    [747](file:///C:/users/peter/github_repos/opty/opty/utils.py:747)                'compilation likely failed. STDERR output from '
    [748](file:///C:/users/peter/github_repos/opty/opty/utils.py:748)                'compilation:\n{}')
--> [749](file:///C:/users/peter/github_repos/opty/opty/utils.py:749)         raise ImportError(msg.format(d['file_prefix'], stderr)) from error
    [750](file:///C:/users/peter/github_repos/opty/opty/utils.py:750) finally:
    [751](file:///C:/users/peter/github_repos/opty/opty/utils.py:751)     module_counter += 1

ImportError: Unable to import the compiled Cython module ufuncify_matrix_1, compilation likely failed. STDERR output from compilation:
c:\Users\Peter\anaconda3\envs\opty-dev\Lib\site-packages\Cython\Compiler\Main.py:381: FutureWarning: Cython directive 'language_level' not set, using '3str' for now (Py3). This has changed from earlier releases! File: C:\Users\Peter\AppData\Local\Temp\tmppwtljn7v.opty_ufuncify_compile\ufuncify_matrix_1.pyx
  tree = Parsing.p_module(s, pxd, full_module_name)
error: command 'C:\\Program Files\\Microsoft Visual Studio\\2022\\Community\\VC\\Tools\\MSVC\\14.38.33130\\bin\\HostX86\\x64\\cl.exe' failed with exit code 2

[moorepants]

That error just tells you that it didn't compile. You have to compile it manually to show what the actual compilation error is.

[tjstienstra]

You have to compile it manually to show what the actual compilation error is.

... >python ufuncify_matrix_1_setup.py build_ext --inplace
ufuncify_matrix_1.c
...\site-packages\numpy\_core\include\numpy\npy_1_7_deprecated_api.h(14) : Warning Msg: Using deprecated NumPy API, disable it with #define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION
ufuncify_matrix_1_c.c
ufuncify_matrix_1_c.c(10): error C2065: 'M_PI': undeclared identifier
error: command '...\\Microsoft Visual Studio\\2022\\Community\\VC\\Tools\\MSVC\\14.40.33807\\bin\\HostX86\\x64\\cl.exe' failed with exit code 2

[moorepants]

I think this issue has already been raised before.

[moorepants]

You opened the same issue on PyDy: pydy/pydy#499

Timo gives the fix there. I suspect it is the same here. It seems that on Windows math constants are not available by default, as they are on linux.

[moorepants]

From https://stackoverflow.com/a/26065595:

Math Constants are not defined in Standard C/C++. To use them, you must first define _USE_MATH_DEFINES and then include cmath or math.h.

[moorepants]

Also see this note in SymPy's C printer: https://github.com/sympy/sympy/blob/d293133e81194adc11177729af91c970f092a6e7/sympy/printing/c.py#L94

[Peter230655]

You opened the same issue on PyDy: pydy/pydy#499

Timo gives the fix there. I suspect it is the same here. It seems that on Windows math constants are not available by default, as they are on linux.

I vaguely remembered that I had raised it before.
I had no intention to make this a big issue again, after all, one fix is easy: replace sm.pi with np.pi.
Main 'strange thing' for me was that sm.pi was not 'accepted', but sm.cos(sm.pi) caused no problems.