Repository: sympy/sympy
Pull request: Raise warning if relation of frames is looped.
<!-- Your title above should be a short description of what
was changed. Do not include the issue number in the title. -->

#### 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. -->
Fixes #21036

#### Brief description of what is fixed or changed
Do bfs to find loop and raise warning.

#### 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.

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 -->
* physics.vector
    * Warning is raised if orientation of frames is looped.
<!-- END RELEASE NOTES -->


Problem statement:
Support orienting adjacent reference frames in arbitrary orders
Suppose you want to establish relative orientation among frames A, B, C, D, and E as such:

```
A ----- B
|
|-------C----D
        |
        |----E
```

A is the root of the tree, B, D, and E are leaves. You do this now with code that looks like:

```python
B.orient(A)
C.orient(A)
D.orient(C)
E.orient(C)
```

This will establish rotation matrices for each connection in the above graph. But a user may, for whatever reason, do something like (call this the alternative use case):

```python
B.orient(A)
A.orient(C)
C.orient(D)
E.orient(C)
```
This currently does not work because with every call of `X.orient()` all adjacent relationships to `X` will be cleared. That is, calling `.orient()` assumes `self`'s orientation relationships should be overwritten. This is sometimes needed, for example if I do:

```python
B.orient(A, ...)
A.orient(B, ...)
```

Anything that was defined in `B.orient(A, ...)` should be fully replaced with the relationships established with `A.orient(B, ...)` because they may be inconsistent with each other. Or if a user tries to do a loop:

```python
B.orient(A)
C.orient(B)
A.orient(C)
```

The last line should raise an error and say something like "Loops in graph not allowed", but what it does is overwrites all relationships to `A` in the last line, effectively undoing the first line.

The alternative use case should work. There is no reason we shouldn't allow construction of the graph in any sequence. Overwriting the relationships to `self` in calls to `orient()` is unnecessary. I think it was implemented like that because it was easier than checking the graph for consistency in a more thorough way.

I think the relationships between points do not have this issue and you can establish them in any order you like. It would be nice if frames also allowed that.

Here is some code that shows how the method `ReferenceFrame._dcm()` wipes relationships of `self`:

```ipython
IPython 7.21.0 -- An enhanced Interactive Python. Type '?' for help.

In [1]: import sympy as sm

In [2]: import sympy.physics.mechanics as me

In [3]: A, B, C = sm.symbols("A, B, C", cls=me.ReferenceFrame)

In [4]: a, b, c = sm.symbols('a, b, c')

In [5]: B.orient(A, 'Axis', (a, A.x))

In [6]: A._dcm_dict
Out[6]: 
{B: Matrix([
 [1,      0,       0],
 [0, cos(a), -sin(a)],
 [0, sin(a),  cos(a)]])}

In [7]: B._dcm_dict
Out[7]: 
{A: Matrix([
 [1,       0,      0],
 [0,  cos(a), sin(a)],
 [0, -sin(a), cos(a)]])}

In [8]: C._dcm_dict
Out[8]: {}

In [9]: B.orient(C, 'Axis', (b, C.x))

In [10]: A._dcm_dict
Out[10]: {}

In [11]: B._dcm_dict
Out[11]: 
{C: Matrix([
 [1,       0,      0],
 [0,  cos(b), sin(b)],
 [0, -sin(b), cos(b)]])}

In [12]: C._dcm_dict
Out[12]: 
{B: Matrix([
 [1,      0,       0],
 [0, cos(b), -sin(b)],
 [0, sin(b),  cos(b)]])}

In [13]: sm.__version__
Out[13]: '1.7.1'
```

Hints from issue discussion:
@angadhn I've opened an issue for this topic that we discussed last night.
A minimum enhancement here would be to explain in the various `.orient*()` documentation that any call to `X.orient*()` will remove any prior relationships to `X` and that to construct a chain or tree of relative oriented reference frames requires that only un-oriented  frames can be added to the ends of the chain/branches.
I'm adding a reference to a previously closed/rejected pull request here. #13824
It is very similar, however, the solution there was to allow loops in a graph.

> Or if a user tries to do a loop:
> 
> ```python
> B.orient(A)
> C.orient(B)
> A.orient(C)
> ```
> 
> The last line should raise an error and say something like "Loops in graph not allowed", but what it does is overwrites all relationships to `A` in the last line, effectively undoing the first line.
> 
> The alternative use case should work. There is no reason we shouldn't allow construction of the graph in any sequence. Overwriting the relationships to `self` in calls to `orient()` is unnecessary. I think it was implemented like that because it was easier than checking the graph for consistency in a more thorough way.
> 
> I think the relationships between points do not have this issue and you can establish them in any order you like. It would be nice if frames also allowed that.


#21271 fixed almost all of the problems mentioned here. Only this part is left to be done.
@moorepants I was wondering if instead of giving error in case of loops in frames, shouldn't we raise `warning`  as we did in velocity.
> @moorepants I was wondering if instead of giving error in case of loops in frames, shouldn't we raise `warning`  as we did in velocity.

@moorepants So it should be a warning or an error?
I guess we should be consistent with what we did with point relationships.
On #21436 I suggested using `topological_sort` for detection of loops. But given the issues described here, is there any reason to not have a helper function to which one passes the orientation information and then the function works out the order -- perhaps through a topological sort -- that orientations should be applied so you don't run into the issues?
I think something like that is needed. We may even need a graph object to manage these things in these modules.

Commit message: Add `Warns` section to all orient functions' docstring

Patch to review:
diff --git a/sympy/physics/vector/frame.py b/sympy/physics/vector/frame.py
index aa9e21e4b561..b3c7a1547ea0 100644
--- a/sympy/physics/vector/frame.py
+++ b/sympy/physics/vector/frame.py
@@ -4,6 +4,8 @@
 from sympy.physics.vector.vector import Vector, _check_vector
 from sympy.utilities.misc import translate
 
+from warnings import warn
+
 __all__ = ['CoordinateSym', 'ReferenceFrame']
 
 
@@ -554,6 +556,25 @@ def _dcm(self, parent, parent_orient):
             self._dcm_dict = self._dlist[0] = {}
         # Reset the _dcm_cache
             self._dcm_cache = {}
+
+        else:
+        #Check for loops and raise warning accordingly.
+            visited = []
+            queue = list(frames)
+            cont = True #Flag to control queue loop.
+            while queue and cont:
+                node = queue.pop(0)
+                if node not in visited:
+                    visited.append(node)
+                    neighbors = node._dcm_dict.keys()
+                    for neighbor in neighbors:
+                        if neighbor == parent:
+                            warn('Loops are defined among the orientation of frames.' + \
+                                ' This is likely not desired and may cause errors in your calculations.')
+                            cont = False
+                            break
+                        queue.append(neighbor)
+
         # Add the dcm relationship to _dcm_dict
         self._dcm_dict.update({parent: parent_orient.T})
         parent._dcm_dict.update({self: parent_orient})
@@ -579,6 +600,12 @@ def orient_axis(self, parent, axis, angle):
         angle : sympifiable
             Angle in radians by which it the frame is to be rotated.
 
+        Warns
+        ======
+
+        UserWarning
+            If frames are looped.
+
         Examples
         ========
 
@@ -657,6 +684,12 @@ def orient_explicit(self, parent, dcm):
             Direction cosine matrix that specifies the relative rotation
             between the two reference frames.
 
+        Warns
+        ======
+
+        UserWarning
+            If frames are looped.
+
         Examples
         ========
 
@@ -761,6 +794,12 @@ def orient_body_fixed(self, parent, angles, rotation_order):
             Tait-Bryan): zxz, xyx, yzy, zyz, xzx, yxy, xyz, yzx, zxy, xzy, zyx,
             and yxz.
 
+        Warns
+        ======
+
+        UserWarning
+            If frames are looped.
+
         Examples
         ========
 
@@ -773,6 +812,7 @@ def orient_body_fixed(self, parent, angles, rotation_order):
         >>> B = ReferenceFrame('B')
         >>> B1 = ReferenceFrame('B1')
         >>> B2 = ReferenceFrame('B2')
+        >>> B3 = ReferenceFrame('B3')
 
         For example, a classic Euler Angle rotation can be done by:
 
@@ -790,8 +830,8 @@ def orient_body_fixed(self, parent, angles, rotation_order):
 
         >>> B1.orient_axis(N, N.x, q1)
         >>> B2.orient_axis(B1, B1.y, q2)
-        >>> B.orient_axis(B2, B2.x, q3)
-        >>> B.dcm(N)
+        >>> B3.orient_axis(B2, B2.x, q3)
+        >>> B3.dcm(N)
         Matrix([
         [        cos(q2),                            sin(q1)*sin(q2),                           -sin(q2)*cos(q1)],
         [sin(q2)*sin(q3), -sin(q1)*sin(q3)*cos(q2) + cos(q1)*cos(q3),  sin(q1)*cos(q3) + sin(q3)*cos(q1)*cos(q2)],
@@ -871,6 +911,12 @@ def orient_space_fixed(self, parent, angles, rotation_order):
             ``'131'``, or the integer ``131``. There are 12 unique valid
             rotation orders.
 
+        Warns
+        ======
+
+        UserWarning
+            If frames are looped.
+
         Examples
         ========
 
@@ -883,6 +929,7 @@ def orient_space_fixed(self, parent, angles, rotation_order):
         >>> B = ReferenceFrame('B')
         >>> B1 = ReferenceFrame('B1')
         >>> B2 = ReferenceFrame('B2')
+        >>> B3 = ReferenceFrame('B3')
 
         >>> B.orient_space_fixed(N, (q1, q2, q3), '312')
         >>> B.dcm(N)
@@ -895,8 +942,8 @@ def orient_space_fixed(self, parent, angles, rotation_order):
 
         >>> B1.orient_axis(N, N.z, q1)
         >>> B2.orient_axis(B1, N.x, q2)
-        >>> B.orient_axis(B2, N.y, q3)
-        >>> B.dcm(N).simplify()
+        >>> B3.orient_axis(B2, N.y, q3)
+        >>> B3.dcm(N).simplify()
         Matrix([
         [ sin(q1)*sin(q2)*sin(q3) + cos(q1)*cos(q3), sin(q1)*cos(q2), sin(q1)*sin(q2)*cos(q3) - sin(q3)*cos(q1)],
         [-sin(q1)*cos(q3) + sin(q2)*sin(q3)*cos(q1), cos(q1)*cos(q2), sin(q1)*sin(q3) + sin(q2)*cos(q1)*cos(q3)],
@@ -992,6 +1039,12 @@ def orient_quaternion(self, parent, numbers):
             The four quaternion scalar numbers as defined above: ``q0``,
             ``q1``, ``q2``, ``q3``.
 
+        Warns
+        ======
+
+        UserWarning
+            If frames are looped.
+
         Examples
         ========
 
@@ -1098,6 +1151,12 @@ def orient(self, parent, rot_type, amounts, rot_order=''):
             ``'123'`` and integer ``123`` are equivalent, for example. Required
             for ``'Body'`` and ``'Space'``.
 
+        Warns
+        ======
+
+        UserWarning
+            If frames are looped.
+
         """
 
         _check_frame(parent)
diff --git a/sympy/physics/vector/tests/test_frame.py b/sympy/physics/vector/tests/test_frame.py
index 61fe2030bbe2..819ba0cb54b2 100644
--- a/sympy/physics/vector/tests/test_frame.py
+++ b/sympy/physics/vector/tests/test_frame.py
@@ -6,6 +6,7 @@
 from sympy.physics.vector.frame import _check_frame
 from sympy.physics.vector.vector import VectorTypeError
 from sympy.testing.pytest import raises
+import warnings
 
 Vector.simp = True
 
@@ -472,6 +473,22 @@ def test_orient_quaternion():
     B.orient_quaternion(A, (0,0,0,0))
     assert B.dcm(A) == Matrix([[0, 0, 0], [0, 0, 0], [0, 0, 0]])
 
+def test_looped_frame_warning():
+    A = ReferenceFrame('A')
+    B = ReferenceFrame('B')
+    C = ReferenceFrame('C')
+
+    a, b, c = symbols('a b c')
+    B.orient_axis(A, A.x, a)
+    C.orient_axis(B, B.x, b)
+
+    with warnings.catch_warnings(record = True) as w:
+        warnings.simplefilter("always")
+        A.orient_axis(C, C.x, c)
+        assert issubclass(w[-1].category, UserWarning)
+        assert 'Loops are defined among the orientation of frames. ' + \
+            'This is likely not desired and may cause errors in your calculations.' in str(w[-1].message)
+
 def test_frame_dict():
     A = ReferenceFrame('A')
     B = ReferenceFrame('B')


Code review comments:
diff_hunk: @@ -579,6 +600,12 @@ def orient_axis(self, parent, axis, angle):
         angle : sympifiable
             Angle in radians by which it the frame is to be rotated.
 
+        Warns
+        ======
+
+        UserWarning
+            If frames are looped.
path: sympy/physics/vector/frame.py
text: @user1:
If the orientation creates a kinematic loop.