Repository: sympy/sympy
Pull request: use only relative error in control_plots test
<!-- Your title above should be a short description of what
was changed. Do not include the issue number in the title. -->
fixes #26728

#### References to other Issues or PRs
<!-- If this pull request fixes an issue, write "Fixes #NNNN" in that exact
format, e.g. "Fixes #1234" (see
https://tinyurl.com/auto-closing for more information). Also, please
write a comment on that issue linking back to this pull request once it is
open. -->


#### Brief description of what is fixed or changed


#### Other comments


#### Release Notes

<!-- Write the release notes for this release below between the BEGIN and END
statements. The basic format is a bulleted list with the name of the subpackage
and the release note for this PR. For example:

* solvers
  * Added a new solver for logarithmic equations.

* functions
  * Fixed a bug with log of integers. Formerly, `log(-x)` incorrectly gave `-log(x)`.

* physics.units
  * Corrected a semantical error in the conversion between volt and statvolt which
    reported the volt as being larger than the statvolt.

or if no release note(s) should be included use:

NO ENTRY

See https://github.com/sympy/sympy/wiki/Writing-Release-Notes for more
information on how to write release notes. The bot will check your release
notes automatically to see if they are formatted correctly. -->

<!-- BEGIN RELEASE NOTES -->
NO ENTRY
<!-- END RELEASE NOTES -->


Problem statement:
Control test fails with numpy 2.0
CC @akshanshbhatt @faze-geek 

This is the second time I have asked to fix these tests (gh-23128).

This test fails with numpy 2.0:
```console
$ pytest sympy/physics/control/tests/test_control_plots.py::test_pole_zero
================================================================================================ test session starts =================================================================================================
platform darwin -- Python 3.12.2, pytest-8.1.1, pluggy-1.4.0
architecture: 64-bit
cache:        yes
ground types: flint (python-flint==0.6.0)

rootdir: /Users/enojb/work/dev/sympy
configfile: pyproject.toml
plugins: instafail-0.5.0, doctestplus-1.2.2.dev5+g5176f32, hypothesis-6.99.9, xdist-3.5.0, timeout-2.3.1, split-0.8.2
collected 1 item

sympy/physics/control/tests/test_control_plots.py F                                                                                                                                                            [100%]

====================================================================================================== FAILURES ======================================================================================================
___________________________________________________________________________________________________ test_pole_zero ___________________________________________________________________________________________________

    def test_pole_zero():
        if not numpy:
            skip("NumPy is required for this test")

        def pz_tester(sys, expected_value):
            z, p = pole_zero_numerical_data(sys)
            z_check = numpy.allclose(z, expected_value[0])
            p_check = numpy.allclose(p, expected_value[1])
            return p_check and z_check

        exp1 = [[], [-0.24999999999999994+1.3919410907075054j, -0.24999999999999994-1.3919410907075054j]]
        exp2 = [[0.0], [-0.25+0.3227486121839514j, -0.25-0.3227486121839514j]]
        exp3 = [[0.0], [-0.5000000000000004+0.8660254037844395j,
            -0.5000000000000004-0.8660254037844395j, 0.9999999999999998+0j]]
        exp4 = [[], [5.0, 0.0, 0.0, 0.0]]
        exp5 = [[-5.645751311064592, -0.5000000000000008, -0.3542486889354093],
            [-0.24999999999999986+1.3919410907075052j,
            -0.24999999999999986-1.3919410907075052j, -0.2499999999999998+0.32274861218395134j,
            -0.2499999999999998-0.32274861218395134j]]
        exp6 = [[], [-1.1641600331447917-3.545808351896439j,
              -0.8358399668552097+2.5458083518964383j]]

>       assert pz_tester(tf1, exp1)
E       assert False
E        +  where False = <function test_pole_zero.<locals>.pz_tester at 0x1148e6c00>(TransferFunction(1, p**2 + 0.5*p + 2, p), [[], [(-0.24999999999999994+1.3919410907075054j), (-0.24999999999999994-1.3919410907075054j)]])

sympy/physics/control/tests/test_control_plots.py:121: AssertionError
                                                                                                   DO *NOT* COMMIT!
============================================================================================== short test summary info ===============================================================================================
FAILED sympy/physics/control/tests/test_control_plots.py::test_pole_zero - assert False
================================================================================================= 1 failed in 0.29s ==================================================================================================
```

Hints from issue discussion:
I'm not able to reproduce this locally. Tests seem to pass with `numpy==2.0.0`.

```
$ python -c "import numpy; print(numpy.version.version)"
2.0.0

$ pytest sympy/physics/control/tests/test_control_plots.py::test_pole_zero
======================== test session starts ========================
platform darwin -- Python 3.12.4, pytest-8.2.2, pluggy-1.5.0
architecture: 64-bit
cache:        yes
ground types: python

rootdir: /Users/akshansh/Documents/GitHub/sympy
configfile: pyproject.toml
plugins: doctestplus-1.2.1, hypothesis-6.103.2, timeout-2.3.1, split-0.8.2, xdist-3.6.1
collected 1 item

sympy/physics/control/tests/test_control_plots.py .           [100%]

========================= 1 passed in 0.16s =========================
```
> I'm not able to reproduce this locally. Tests seem to pass with `numpy==2.0.0`.

I can reproduce this only on MacOS. On Linux it passes.
On MacOS what is returned is:
```
(Pdb) p pole_zero_numerical_data(tf1)
(array([], dtype=float64), array([-0.25-1.39194109j, -0.25+1.39194109j]))
```
The test expects
```
        exp1 = [[], [-0.24999999999999994+1.3919410907075054j, -0.24999999999999994-1.3919410907075054j]]
```
The difference is the ordering of the roots.
The underlying call is to np.roots:
```python
In [1]: np.roots([1. +0.j, 0.5+0.j, 2. +0.j]) # numpy 1.26.4 on MacOS
Out[1]: array([-0.25+1.39194109j, -0.25-1.39194109j])
```
```python
In [1]: np.roots([1. +0.j, 0.5+0.j, 2. +0.j]) # numpy 2.0.0 on MacOS
Out[1]: array([-0.25-1.39194109j, -0.25+1.39194109j])
```
The documentation for np.roots does not say anything about the ordering of the returned roots. It says that they come from computing the eigenvalues of the companion matrix. The eigvals function does say that the roots are not necessarily ordered.

The docs do say:
```
    .. note::
       This forms part of the old polynomial API. Since version 1.4, the
       new polynomial API defined in `numpy.polynomial` is preferred.
       A summary of the differences can be found in the
       :doc:`transition guide </reference/routines.polynomials>`.
```
Really though we should be using sympy routines to find the roots because np.roots can be inaccurate for some roots.
This is the guaranteed accurate way to find numerical approximations of the roots for any polynomial with rational or float coefficients using SymPy's routines:
```python
In [20]: p = x**2 + x/2 + 2

In [21]: [r.evalf() for r in all_roots(p)]
Out[21]: [-0.25 - 1.39194109070751⋅ⅈ, -0.25 + 1.39194109070751⋅ⅈ]
```
This method uses exact arithmetic to isolate the roots, to distinguish real and non-real roots, and to identify complex conjugate pairs. The roots are returned by `all_roots` in a guaranteed order that is real roots followed by non-real roots in pairs ordered first by real parts and then by the magnitudes of the imaginary parts:
```python
In [22]: all_roots((x-1)*(x-I)*(x+I)*(x-2*I)*(x+2*I))
Out[22]: [1, -ⅈ, ⅈ, -2⋅ⅈ, 2⋅ⅈ]
```
The `all_roots` function is not limited by degree because it uses `RootOf`. The call to `evalf` from the exact representation of the roots is guaranteed to give a result that is accurate to the requested precision.

NumPy's `np.roots` function cannot provide similar guarantees because it uses fixed precision floating point.
NumPy would ideally be eliminated from `control_plots.py`. This code gets merged without being stricter about eating our own dog food for plotting.
I understand that using numpy is not the most elegant approach to get polynomial root data. Still, I don't think that sympy's built-in polynomial root calculator is fast enough to use it instead. Sure, sympy gives more consistent and precise results, but we don't require low-tolerance data points for plotting purposes. See the time difference using both of these methods:

```python
In [2]: p = x**8 + 9*x**5 + 11*x + 7  # Arbitrary 8-deg polynomial

In [3]: %time [r.evalf() for r in all_roots(p)]
CPU times: user 4.1 s, sys: 5.52 ms, total: 4.11 s
Wall time: 4.11 s
Out[3]:
[-2.10954725209519,
 -0.582637006756454,
 -0.55581013748193 - 0.805706304516959*I,
 -0.55581013748193 + 0.805706304516959*I,
 0.912751454705842 - 0.742665591979299*I,
 0.912751454705842 + 0.742665591979299*I,
 0.98915081220191 - 1.82059670736243*I,
 0.98915081220191 + 1.82059670736243*I]

In [4]: %time np.roots(np.array(Poly(p, x).all_coeffs(), dtype=np.complex128))
CPU times: user 2.41 ms, sys: 5.06 ms, total: 7.47 ms
Wall time: 11.9 ms
Out[4]:
array([ 0.98915081+1.82059671e+00j,  0.98915081-1.82059671e+00j,
       -2.10954725+3.95682719e-16j,  0.91275145+7.42665592e-01j,
        0.91275145-7.42665592e-01j, -0.55581014+8.05706305e-01j,
       -0.55581014-8.05706305e-01j, -0.58263701-2.23618104e-25j])
```

The precision vs calculation-time tradeoff is quite skewed. This is the reason we decided to use numpy in the first place. Sometime in the future, if root-locus plot is added to this module, it will depend on the pz numerical value function for root calculations, which will be called a hundred, perhaps even a thousand times (for sampling loci points). It will only make this time bottleneck worse. I would love to hear if you have any suggestions for solving this.
We can make the SymPy routines faster in various ways including by using better algorithms and/or using python-flint:
```python
In [3]: from flint import *

In [4]: x = fmpz_poly([0, 1])

In [5]: p = x**8 + 9*x**5 + 11*x + 7

In [6]: p
Out[6]: x^8 + 9*x^5 + 11*x + 7

In [7]: %time p.complex_roots()
CPU times: user 6.61 ms, sys: 0 ns, total: 6.61 ms
Wall time: 76.5 ms
Out[7]: 
[([-2.10954725209519 +/- 6.57e-16], 1),
 ([-0.582637006756454 +/- 3.76e-16], 1),
 ([0.912751454705842 +/- 4.04e-16] + [0.742665591979299 +/- 1.98e-16]j, 1),
 ([0.912751454705842 +/- 4.04e-16] + [-0.742665591979299 +/- 1.98e-16]j, 1),
 ([-0.555810137481930 +/- 4.43e-16] + [0.805706304516960 +/- 4.93e-16]j, 1),
 ([-0.555810137481930 +/- 4.43e-16] + [-0.805706304516960 +/- 4.93e-16]j, 1),
 ([0.989150812201910 +/- 4.79e-16] + [1.82059670736243 +/- 1.14e-15]j, 1),
 ([0.989150812201910 +/- 4.79e-16] + [-1.82059670736243 +/- 1.14e-15]j, 1)]
```

> The precision vs calculation-time tradeoff is quite skewed. This is the reason we decided to use numpy in the first place

Why should anyone use SymPy's control module at all?

There are other Python control modules. If someone wants to use numpy to do a root locus plot then they can use the `control` module:
https://python-control.readthedocs.io/en/latest/generated/control.root_locus.html

What purpose is served by SymPy reimplementing things that already exist in other libraries?

Providing more accurate calculations is one possible reason for SymPy to have a control module. Duplicating the features that already exist in other Python libraries is not.

If SymPy's control module plotting functionality is not any better than the other control libraries then I don't think any plotting functionality should be added to SymPy's control module apart from just calling into the other libraries.
> ```python
> In [7]: %time p.complex_roots()
> CPU times: user 6.61 ms, sys: 0 ns, total: 6.61 ms
> Wall time: 76.5 ms
> ```

Note actually that this time was inflated by initialising the printer. Here is the actual time:
```python
In [1]: from flint import *

In [2]: x = fmpz_poly([0, 1])

In [3]: p = x**8 + 9*x**5 + 11*x + 7

In [4]: %time r = p.complex_roots()
CPU times: user 2.19 ms, sys: 0 ns, total: 2.19 ms
Wall time: 2.21 ms
```

Commit message: Update test_control_plots.py

Patch to review:
diff --git a/sympy/physics/control/control_plots.py b/sympy/physics/control/control_plots.py
index 3742de329e61..7f1fc5179ee6 100644
--- a/sympy/physics/control/control_plots.py
+++ b/sympy/physics/control/control_plots.py
@@ -2,12 +2,14 @@
 from sympy.functions.elementary.exponential import (exp, log)
 from sympy.polys.partfrac import apart
 from sympy.core.symbol import Dummy
+from sympy.core.sympify import _sympify
 from sympy.external import import_module
 from sympy.functions import arg, Abs
 from sympy.integrals.laplace import _fast_inverse_laplace
 from sympy.physics.control.lti import SISOLinearTimeInvariant
 from sympy.plotting.series import LineOver1DRangeSeries
 from sympy.polys.polytools import Poly
+from sympy.polys.polyutils import _nsort
 from sympy.printing.latex import latex
 
 __all__ = ['pole_zero_numerical_data', 'pole_zero_plot',
@@ -88,8 +90,8 @@ def pole_zero_numerical_data(system):
     >>> from sympy.physics.control.lti import TransferFunction
     >>> from sympy.physics.control.control_plots import pole_zero_numerical_data
     >>> tf1 = TransferFunction(s**2 + 1, s**4 + 4*s**3 + 6*s**2 + 5*s + 2, s)
-    >>> pole_zero_numerical_data(tf1)   # doctest: +SKIP
-    ([-0.+1.j  0.-1.j], [-2. +0.j        -0.5+0.8660254j -0.5-0.8660254j -1. +0.j       ])
+    >>> pole_zero_numerical_data(tf1)
+    ([-1j, 1j], [-2.0, (-1.5-0.8660254j), (-0.5+0.8660254j)])
 
     See Also
     ========
@@ -109,7 +111,12 @@ def pole_zero_numerical_data(system):
     zeros = np.roots(num_poly)
     poles = np.roots(den_poly)
 
-    return zeros, poles
+    # make ordering canonical
+    def _sort(l):
+        return [float(i) if i.is_real else complex(i) for i in
+                _nsort([_sympify(i) for i in l])]
+
+    return _sort(zeros), _sort(poles)
 
 
 def pole_zero_plot(system, pole_color='blue', pole_markersize=10,
diff --git a/sympy/physics/control/tests/test_control_plots.py b/sympy/physics/control/tests/test_control_plots.py
index 673fcee6cfdb..eeb73e6d3702 100644
--- a/sympy/physics/control/tests/test_control_plots.py
+++ b/sympy/physics/control/tests/test_control_plots.py
@@ -1,5 +1,6 @@
 from math import isclose
-from sympy.core.numbers import I
+from sympy.core.numbers import I, all_close
+from sympy.core.sympify import sympify
 from sympy.core.symbol import Dummy
 from sympy.functions.elementary.complexes import (Abs, arg)
 from sympy.functions.elementary.exponential import log
@@ -101,20 +102,27 @@ def test_pole_zero():
         skip("NumPy is required for this test")
 
     def pz_tester(sys, expected_value):
-        z, p = pole_zero_numerical_data(sys)
-        z_check = numpy.allclose(z, expected_value[0])
-        p_check = numpy.allclose(p, expected_value[1])
+        z, p = [list(i) for i in pole_zero_numerical_data(sys)]
+        def check(a, b):
+            if isinstance(a, (list, tuple)):
+                return all(check(i, j) for i,j in zip(a, b))
+            a = sympify(a).n(chop=1e-10)
+            if not b:
+                return all_close(a, b, atol=1e-10, rtol=0)
+            return all_close(a, b, atol=0, rtol=1e-6)
+        z_check = check(z, expected_value[0])
+        p_check = check(p, expected_value[1])
         return p_check and z_check
 
-    exp1 = [[], [-0.24999999999999994+1.3919410907075054j, -0.24999999999999994-1.3919410907075054j]]
-    exp2 = [[0.0], [-0.25+0.3227486121839514j, -0.25-0.3227486121839514j]]
-    exp3 = [[0.0], [-0.5000000000000004+0.8660254037844395j,
-        -0.5000000000000004-0.8660254037844395j, 0.9999999999999998+0j]]
-    exp4 = [[], [5.0, 0.0, 0.0, 0.0]]
+    exp1 = [[], [-0.24999999999999994-1.3919410907075054j, -0.24999999999999994+1.3919410907075054j]]
+    exp2 = [[0.0], [-0.25-0.3227486121839514j, -0.25+0.3227486121839514j]]
+    exp3 = [[0.0], [0.9999999999999998+0j, -0.5000000000000004-0.8660254037844395j,
+        -0.5000000000000004+0.8660254037844395j]]
+    exp4 = [[], [0.0, 0.0, 0.0, 5.0]]
     exp5 = [[-5.645751311064592, -0.5000000000000008, -0.3542486889354093],
-        [-0.24999999999999986+1.3919410907075052j,
-        -0.24999999999999986-1.3919410907075052j, -0.2499999999999998+0.32274861218395134j,
-        -0.2499999999999998-0.32274861218395134j]]
+        [-0.24999999999999986-1.3919410907075052j,
+        -0.24999999999999986-0.32274861218395134j, -0.2499999999999998+0.32274861218395134j,
+        -0.2499999999999998+1.3919410907075052j]]
     exp6 = [[], [-1.1641600331447917-3.545808351896439j,
           -0.8358399668552097+2.5458083518964383j]]
 


Code review comments:
diff_hunk: @@ -88,8 +90,8 @@ def pole_zero_numerical_data(system):
     >>> from sympy.physics.control.lti import TransferFunction
     >>> from sympy.physics.control.control_plots import pole_zero_numerical_data
     >>> tf1 = TransferFunction(s**2 + 1, s**4 + 4*s**3 + 6*s**2 + 5*s + 2, s)
-    >>> pole_zero_numerical_data(tf1)   # doctest: +SKIP
-    ([-0.+1.j  0.-1.j], [-2. +0.j        -0.5+0.8660254j -0.5-0.8660254j -1. +0.j       ])
+    >>> pole_zero_numerical_data(tf1)
path: sympy/physics/control/control_plots.py
text: @author:
```suggestion
    >>> pole_zero_numerical_data(tf1)  # doctest: +SKIP
```