Repository: scipy/scipy
Pull request: BUG: Fix error with 180 degree rotation in Rotation.align_vectors() with an infinite weight
#### Reference issue
Closes https://github.com/scipy/scipy/issues/20555

#### What does this implement/fix?
There are two related issues with exact rotations near 180 degrees:

* A perfect 180 deg rotation results in r = cross = [0, 0, 0], and is thus cast to a 0 degree rotation (critical defect)
* Numerical instabilities near a 180 deg rotation (non-critical defect)

This fixes both. For the numerical instability part, near 180 degrees for the test case I'm seeing a pointing error of 2e-6 drop to 5e-9 (and now matches the non-exact code path to 1e-11). Tried to make the code path as compact as possible, but improvements welcome!

Problem statement:
BUG: spatial: error in `Rotation.align_vectors()` with an infinite weight
### Describe your issue.

I found a degenerate case when aligning vectors exactly with `Rotation.align_vectors` and an infinite weight. 

### Reproducing Code Example

```python
import numpy as np
from scipy.spatial.transform import Rotation

# Wrong:
Rotation.align_vectors([[0, 1, 0], [-1, 0, 0]], [[0, -1, 0], [-1, 0, 0]], weights=[np.inf, 1])[0].as_rotvec()
# array([0., 0., 0.])

# Right:
Rotation.align_vectors([[0, 1, 0], [-1, 0, 0]], [[0, -1, 0], [-1, 0, 0]], weights=[1e6, 1])[0].as_rotvec()
# array([3.14159265, 0., 0.])
```

### SciPy/NumPy/Python version and system information

```shell
import sys, scipy, numpy; print(scipy.__version__, numpy.__version__, sys.version_info); scipy.show_config()
1.13.0 1.26.3 sys.version_info(major=3, minor=10, micro=12, releaselevel='final', serial=0)
Build Dependencies:
  blas:
    detection method: pkgconfig
    found: true
    include directory: /usr/local/include
    lib directory: /usr/local/lib
    name: openblas
    openblas configuration: USE_64BITINT=0 DYNAMIC_ARCH=1 DYNAMIC_OLDER= NO_CBLAS=
      NO_LAPACK= NO_LAPACKE= NO_AFFINITY=1 USE_OPENMP= ZEN MAX_THREADS=64
    pc file directory: /usr/local/lib/pkgconfig
    version: 0.3.26.dev
  lapack:
    detection method: pkgconfig
    found: true
    include directory: /usr/local/include
    lib directory: /usr/local/lib
    name: openblas
    openblas configuration: USE_64BITINT=0 DYNAMIC_ARCH=1 DYNAMIC_OLDER= NO_CBLAS=
      NO_LAPACK= NO_LAPACKE= NO_AFFINITY=1 USE_OPENMP= ZEN MAX_THREADS=64
    pc file directory: /usr/local/lib/pkgconfig
    version: 0.3.26.dev
  pybind11:
    detection method: config-tool
    include directory: unknown
    name: pybind11
    version: 2.12.0
Compilers:
  c:
    commands: cc
    linker: ld.bfd
    name: gcc
    version: 10.2.1
  c++:
    commands: c++
    linker: ld.bfd
    name: gcc
    version: 10.2.1
  cython:
    commands: cython
    linker: cython
    name: cython
    version: 3.0.10
  fortran:
    commands: gfortran
    linker: ld.bfd
    name: gcc
    version: 10.2.1
  pythran:
    include directory: ../../tmp/pip-build-env-0blqy1or/overlay/lib/python3.10/site-packages/pythran
    version: 0.15.0
Machine Information:
  build:
    cpu: x86_64
    endian: little
    family: x86_64
    system: linux
  cross-compiled: false
  host:
    cpu: x86_64
    endian: little
    family: x86_64
    system: linux
Python Information:
  path: /opt/python/cp310-cp310/bin/python
  version: '3.10'
```

Hints from issue discussion:
I'm seeing the issue here, working on a fix. Problem is here:
https://github.com/scipy/scipy/blob/d55cb95282b5cfab381e373d3c127630c7e915a0/scipy/spatial/transform/_rotation.pyx#L3483-L3492

There are two related problems:
* A perfect 180 deg rotation results in `r = cross = [0, 0, 0]` (this bug, an error)
* Potential numerical instabilities near a 180 deg rotation (not exactly this bug but will fix and add tests)

Commit message: Tests for exact near 180 deg rotations

Patch to review:
diff --git a/scipy/spatial/transform/_rotation.pyx b/scipy/spatial/transform/_rotation.pyx
index f18a24270277..a1e29747ab1c 100644
--- a/scipy/spatial/transform/_rotation.pyx
+++ b/scipy/spatial/transform/_rotation.pyx
@@ -1243,7 +1243,7 @@ cdef class Rotation:
         for ind in range(num_rotations):
             angle = _norm3(crotvec[ind, :])
 
-            if angle <= 1e-3:  # small angle Taylor series expansion
+            if angle <= 1e-3:  # small angle order 5 Taylor series expansion
                 angle2 = angle * angle
                 scale = 0.5 - angle2 / 48 + angle2 * angle2 / 3840
             else:  # large angle
@@ -1916,7 +1916,7 @@ cdef class Rotation:
 
             angle = 2 * atan2(_norm3(quat), quat[3])
 
-            if angle <= 1e-3:  # small angle Taylor series expansion
+            if angle <= 1e-3:  # small angle order 5 Taylor series expansion
                 angle2 = angle * angle
                 scale = 2 + angle2 / 12 + 7 * angle2 * angle2 / 2880
             else:  # large angle
@@ -3483,13 +3483,32 @@ cdef class Rotation:
             # We first find the minimum angle rotation between the primary
             # vectors.
             cross = np.cross(b_pri[0], a_pri[0])
-            theta = atan2(_norm3(cross), np.dot(a_pri[0], b_pri[0]))
-            if theta < 1e-3:  # small angle Taylor series approximation
+            cross_norm = _norm3(cross)
+            theta = atan2(cross_norm, np.dot(a_pri[0], b_pri[0]))
+            tolerance = 1e-3  # tolerance for small angle approximation (rad)
+            R_flip = cls.identity()
+            if abs(np.pi - theta) < tolerance:
+                # At 180 degrees, cross = [0, 0, 0] so we need to flip the
+                # vectors (along an arbitrary orthogonal axis of rotation)
+                if cross_norm == 0:
+                    vec_non_parallel = np.roll(a_pri, 1) + np.array([1, 0, 0])
+                    vec_orthogonal = np.cross(a_pri, vec_non_parallel)
+                    r_flip = vec_orthogonal / _norm3(vec_orthogonal) * np.pi
+                    R_flip = cls.from_rotvec(r_flip)
+
+                # Near 180 degrees, the small angle appoximation of x/sin(x)
+                # diverges, so for numerical stability we flip and then rotate
+                # back by pi - theta
+                else:
+                    R_flip = cls.from_rotvec(cross / _norm3(cross) * np.pi)
+                    theta = np.pi - theta
+                    cross = -cross
+            if theta < 1e-3:  # small angle order 5 Taylor series approximation
                 theta2 = theta * theta
                 r = cross * (1 + theta2 / 6 + theta2 * theta2 * 7 / 360)
             else:
                 r = cross * theta / np.sin(theta)
-            R_pri = cls.from_rotvec(r)
+            R_pri = cls.from_rotvec(r) * R_flip
 
             if N == 1:
                 # No secondary vectors, so we are done
diff --git a/scipy/spatial/transform/tests/test_rotation.py b/scipy/spatial/transform/tests/test_rotation.py
index bf208f34437b..bc9af0363ae0 100644
--- a/scipy/spatial/transform/tests/test_rotation.py
+++ b/scipy/spatial/transform/tests/test_rotation.py
@@ -1467,6 +1467,24 @@ def test_align_vectors_parallel():
     assert_allclose(R.apply(b[0]), a[0], atol=atol)
 
 
+def test_align_vectors_antiparallel():
+    # Test exact 180 deg rotation
+    atol = 1e-12
+    a = [[0, 1, 0], [-1, 0, 0]]
+    b = [[0, -1, 0], [-1, 0, 0]]
+    R, _ = Rotation.align_vectors(a, b, weights=[np.inf, 1])
+    assert np.isclose(R.magnitude(), np.pi, atol=atol)
+
+    # Test exact rotations near 180 deg
+    as_to_test = [a,
+                  [[0, 1,  1e-4], [-1, 0, 0]],
+                  [[0, 1, -1e-4], [-1, 0, 0]]]
+    for a in as_to_test:
+        R, _ = Rotation.align_vectors(a, b, weights=[np.inf, 1])
+        R2, _ = Rotation.align_vectors(a, b, weights=[1e10, 1])
+        assert R.approx_equal(R2, atol=atol)
+
+
 def test_align_vectors_primary_only():
     atol = 1e-12
     mats_a = Rotation.random(100, random_state=0).as_matrix()


Code review comments:
diff_hunk: @@ -1467,6 +1467,24 @@ def test_align_vectors_parallel():
     assert_allclose(R.apply(b[0]), a[0], atol=atol)
 
 
+def test_align_vectors_antiparallel():
+    # Test exact 180 deg rotation
+    atol = 1e-12
+    a = [[0, 1, 0], [-1, 0, 0]]
+    b = [[0, -1, 0], [-1, 0, 0]]
+    R, _ = Rotation.align_vectors(a, b, weights=[np.inf, 1])
+    assert np.isclose(R.magnitude(), np.pi, atol=atol)
path: scipy/spatial/transform/tests/test_rotation.py
text: @user1:
Just a surface-level comment for now, but could you use the built-in NumPy testing functions for making these kinds of assertions? The feedback is much better, even for scalars:

```python
>>> import numpy as np
>>> from numpy.testing import assert_allclose
>>> assert np.isclose(3.1, np.pi)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
AssertionError
>>> assert_allclose(3.1, np.pi)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/home/treddy/python_venvs/py_311_scipy_dev/lib/python3.11/site-packages/numpy/testing/_private/utils.py", line 1684, in assert_allclose
    assert_array_compare(compare, actual, desired, err_msg=str(err_msg),
  File "/home/treddy/github_projects/spack/opt/spack/linux-ubuntu22.04-skylake/gcc-11.3.0/python-3.11.2-4qncg2nqev5evxfmamde3e6rnb34b4ls/lib/python3.11/contextlib.py", line 81, in inner
    return func(*args, **kwds)
           ^^^^^^^^^^^^^^^^^^^
  File "/home/treddy/python_venvs/py_311_scipy_dev/lib/python3.11/site-packages/numpy/testing/_private/utils.py", line 885, in assert_array_compare
    raise AssertionError(msg)
AssertionError: 
Not equal to tolerance rtol=1e-07, atol=0

Mismatched elements: 1 / 1 (100%)
Max absolute difference among violations: 0.04159265
Max relative difference among violations: 0.01323935
 ACTUAL: array(3.1)
 DESIRED: array(3.141593)
```