Repository: sympy/sympy
Pull request: fix(diffgeom): fix indirect transformation
<!-- 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 #20742

#### Brief description of what is fixed or changed

Allow indirect transformation
Change signature of relation in order to allow indirect transformation
#### Other comments

Before #21464, transformation relation for `CoordSystem` used `Lambda` e.g. `Lambda((x, y), (y, -x))`.
In #21464, the signature was changed to use only the expression, e.g. `(y, -x)`.

However, to allow indirect transformation it is crucial to store the original symbol and transformed expression altogether. Therefore, the signature is now `((x, y), (y, -x))`. Old signature with `Lambda` is still supported, so it won't break backwards compatibility.

#### 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 -->
* diffgeom
    * Bug in indirect transformation between coordinate systems is fixed.
    * Transform relation between coordinate systems now consists of two tuples instead of `Lambda`.
<!-- END RELEASE NOTES -->


Problem statement:
diffgeom.CoordSystem.transform with indirect relation raises KeyError
Below code is a minimal example for an indirect transform via a temporary intermediate coordinate system. `parabolic2D.transform(poloidal)` raises a KeyError.

MWE:
```
"""
2D manifold coordinate systems
"""
import sympy
import sympy.diffgeom

CARTESIAN2D = "cartesian2D"
POLOIDAL = "poloidal"
PARABOLIC2D = "parabolic2D"

m = sympy.diffgeom.Manifold("M", 2)
p = sympy.diffgeom.Patch("P", m)
RELATIONS = {}

# #########
# cartesian
# #########
x, y = sympy.symbols("x y")

# #########
# poloidal
# #########
r, theta = sympy.symbols("r theta", nonnegative=True)
RELATIONS[(CARTESIAN2D, POLOIDAL)] = sympy.Lambda(
    (x, y),
    sympy.Matrix(
        [
            sympy.sqrt(x ** 2 + y ** 2),
            sympy.atan2(y, x)
        ]
    )
)
RELATIONS[(POLOIDAL, CARTESIAN2D)] = sympy.Lambda(
    (r, theta),
    sympy.Matrix(
        [
            r * sympy.cos(theta),
            r * sympy.sin(theta)
        ]
    )
)

# #########
# parabolic
# #########
sigma, tau = sympy.symbols("sigma tau")
RELATIONS[(PARABOLIC2D, CARTESIAN2D)] = sympy.Lambda(
    (sigma, tau),
    sympy.Matrix(
        [
            sigma * tau,
            1 / 2 * (tau**2 - sigma**2)
        ]
    )
)

cartesian2D = sympy.diffgeom.CoordSystem(CARTESIAN2D, p, [x, y], RELATIONS)
poloidal = sympy.diffgeom.CoordSystem(POLOIDAL, p, [r, theta], RELATIONS)
parabolic2D = sympy.diffgeom.CoordSystem(PARABOLIC2D, p, [sigma, tau], RELATIONS)


if __name__ == "__main__":
    print(parabolic2D.transform(poloidal))  # raises a KeyError
    print(poloidal.transform(parabolic2D))  # raises a KeyError
```

This raises a KeyError.

> Traceback (most recent call last):
>   File "/opt/anaconda/envs/py38/lib/python3.8/pdb.py", line 1703, in main
>     pdb._runscript(mainpyfile)
>   File "/opt/anaconda/envs/py38/lib/python3.8/pdb.py", line 1572, in _runscript
>     self.run(statement)
>   File "/opt/anaconda/envs/py38/lib/python3.8/bdb.py", line 580, in run
>     exec(cmd, globals, locals)
>   File "<string>", line 1, in <module>
>   File "/home/IPP-HGW/dboe/git/tfields/tfields/bases/manifold_2.py", line 1, in <module>
>     """
>   File "/opt/anaconda/envs/py38/lib/python3.8/site-packages/sympy/diffgeom/diffgeom.py", line 480, in transform
>     transf = self.transformation(sys)
>   File "/opt/anaconda/envs/py38/lib/python3.8/site-packages/sympy/diffgeom/diffgeom.py", line 354, in transformation
>     return self._indirect_transformation(self, sys)
>   File "/opt/anaconda/envs/py38/lib/python3.8/site-packages/sympy/core/cache.py", line 72, in wrapper
>     retval = cfunc(*args, **kwargs)
>   File "/opt/anaconda/envs/py38/lib/python3.8/site-packages/sympy/diffgeom/diffgeom.py", line 376, in _indirect_transformation
>     path = cls._dijkstra(sys1, sys2)
>   File "/opt/anaconda/envs/py38/lib/python3.8/site-packages/sympy/diffgeom/diffgeom.py", line 414, in _dijkstra
>     visit(sys1)
>   File "/opt/anaconda/envs/py38/lib/python3.8/site-packages/sympy/diffgeom/diffgeom.py", line 406, in visit
>     path_dict[sys][2] = 1
> KeyError: parabolic2D
> 

I think I found the reson already: In dijkstra routing the comparison is done between a CoordSystem and sympy.Str
DEBUGGING:

```
Uncaught exception. Entering post mortem debugging
Running 'cont' or 'step' will restart the program
> /opt/anaconda/envs/py38/lib/python3.8/site-packages/sympy/diffgeom/diffgeom.py(406)visit()
-> path_dict[sys][2] = 1
(Pdb) path_dict
{cartesian2D: [0, [], 0], poloidal: [0, [], 0], parabolic2D: [0, [], 0]}
(Pdb) sys
parabolic2D
(Pdb) hash(sys)
-2150956724454717966
(Pdb) [hash(k) for k in path_dict]
[6233656281303402583, 5480353473597806494, -1340528192013030397]
(Pdb) type(sys)
<class 'sympy.diffgeom.diffgeom.CoordSystem'>
(Pdb) [type(k) for k in path_dict]
[<class 'sympy.core.symbol.Str'>, <class 'sympy.core.symbol.Str'>, <class 'sympy.core.symbol.Str'>]
```

Hints from issue discussion:
I started looking into this, I noticed another bug in addition to the type mismatch in the Dijkstra algorithm.  Your example only gives the forward transformation from parabolic2D to cartesian2D, so in order to solve `poloidal.transform(parabolic2D)` the `transform `method will need to calculate the inverse transformation from cartesian2D to parabolic2D, which will raise an exception, see https://github.com/sympy/sympy/issues/21356

Commit message: fix(diffgeom): fix indirect transformation

Allow indirect transformation
Change signature of relation in order to allow indirect transformation

Patch to review:
diff --git a/sympy/diffgeom/diffgeom.py b/sympy/diffgeom/diffgeom.py
index ba82beebbab6..1e7b04be4c95 100644
--- a/sympy/diffgeom/diffgeom.py
+++ b/sympy/diffgeom/diffgeom.py
@@ -193,8 +193,9 @@ class CoordSystem(Basic):
 
     relations : dict, optional
         Key is a tuple of two strings, who are the names of the systems where
-        the coordinates transform from and transform to. Value is a tuple of
-        transformed coordinates.
+        the coordinates transform from and transform to.
+        Value is a tuple of the symbols before transformation and a tuple of
+        the expressions after transformation.
 
     Examples
     ========
@@ -209,8 +210,8 @@ class CoordSystem(Basic):
     >>> x, y = symbols('x y', real=True)
     >>> r, theta = symbols('r theta', nonnegative=True)
     >>> relation_dict = {
-    ... ('Car2D', 'Pol'): (sqrt(x**2 + y**2), atan2(y, x)),
-    ... ('Pol', 'Car2D'): (r*cos(theta), r*sin(theta))
+    ... ('Car2D', 'Pol'): [(x, y), (sqrt(x**2 + y**2), atan2(y, x))],
+    ... ('Pol', 'Car2D'): [(r, theta), (r*cos(theta), r*sin(theta))]
     ... }
     >>> Car2D = CoordSystem('Car2D', p, (x, y), relation_dict)
     >>> Pol = CoordSystem('Pol', p, (r, theta), relation_dict)
@@ -314,8 +315,12 @@ def __new__(cls, name, patch, symbols=None, relations={}, **kwargs):
             if not isinstance(s2, Str):
                 s2 = Str(s2)
             key = Tuple(s1, s2)
+
+            # Old version used Lambda as a value.
             if isinstance(v, Lambda):
-                v = tuple(v(*symbols))
+                v = (tuple(v.signature), tuple(v.expr))
+            else:
+                v = (tuple(v[0]), tuple(v[1]))
             rel_temp[key] = v
         relations = Dict(rel_temp)
 
@@ -396,7 +401,7 @@ def transformation(self, sys):
         if self == sys:
             expr = Matrix(self.symbols)
         elif key in self.relations:
-            expr = Matrix(self.relations[key])
+            expr = Matrix(self.relations[key][1])
         elif key[::-1] in self.relations:
             expr = Matrix(self._inverse_transformation(sys, self))
         else:
@@ -404,44 +409,58 @@ def transformation(self, sys):
         return Lambda(signature, expr)
 
     @staticmethod
-    def _inverse_transformation(sys1, sys2):
-        # Find the transformation relation from sys2 to sys1
-        forward_transform_expressions = sys1.transform(sys2)
-
-        inv_results = solve(
-            [t[0] - t[1] for t in zip(sys2.symbols, forward_transform_expressions)],
-            list(sys1.symbols), dict=True)
-        if len(inv_results) == 0:
-            raise NotImplementedError(
-                "Cannot solve inverse of transformation from {} to {}".format(sys1, sys2))
-        elif len(inv_results) > 1:
-            raise ValueError(
-                "Obtained multiple results for inverse of transformation from {} to {}".format(sys1, sys2)
-            )
+    def _solve_inverse(sym1, sym2, exprs, sys1_name, sys2_name):
+        ret = solve(
+            [t[0] - t[1] for t in zip(sym2, exprs)],
+            list(sym1), dict=True)
+
+        if len(ret) == 0:
+            temp = "Cannot solve inverse relation from {} to {}."
+            raise NotImplementedError(temp.format(sys1_name, sys2_name))
+        elif len(ret) > 1:
+            temp = "Obtained multiple inverse relation from {} to {}."
+            raise ValueError(temp.format(sys1_name, sys2_name))
+
+        return ret[0]
 
-        inv_results = inv_results[0]
+    @classmethod
+    def _inverse_transformation(cls, sys1, sys2):
+        # Find the transformation relation from sys2 to sys1
+        forward = sys1.transform(sys2)
+        inv_results = cls._solve_inverse(sys1.symbols, sys2.symbols, forward,
+                                         sys1.name, sys2.name)
         signature = tuple(sys1.symbols)
         return [inv_results[s] for s in signature]
 
     @classmethod
     @cacheit
     def _indirect_transformation(cls, sys1, sys2):
-        # Find the transformation relation between two indirectly connected coordinate systems
+        # Find the transformation relation betweentwo indirectly connected
+        # coordinate systems
+        rel = sys1.relations
         path = cls._dijkstra(sys1, sys2)
-        Lambdas = []
-        for i in range(len(path) - 1):
-            s1, s2 = path[i], path[i + 1]
-            Lambdas.append(s1.transformation(s2))
-        syms = Lambdas[-1].signature
-        expr = syms
-        for l in reversed(Lambdas):
-            expr = l(*expr)
-        return Lambda(syms, expr)
+
+        transforms = []
+        for s1, s2 in zip(path, path[1:]):
+            if (s1, s2) in rel:
+                transforms.append(rel[(s1, s2)])
+            else:
+                sym2, inv_exprs = rel[(s2, s1)]
+                sym1 = tuple(Dummy() for i in sym2)
+                ret = cls._solve_inverse(sym2, sym1, inv_exprs, s2, s1)
+                ret = tuple(ret[s] for s in sym2)
+                transforms.append((sym1, ret))
+        syms = sys1.args[2]
+        exprs = syms
+        for newsyms, newexprs in transforms:
+            exprs = tuple(e.subs(zip(newsyms, exprs)) for e in newexprs)
+        return exprs
 
     @staticmethod
     def _dijkstra(sys1, sys2):
         # Use Dijkstra algorithm to find the shortest path between two indirectly-connected
         # coordinate systems
+        # return value is the list of the names of the systems.
         relations = sys1.relations
         graph = {}
         for s1, s2 in relations.keys():
@@ -465,7 +484,7 @@ def visit(sys):
                     path_dict[newsys][1] = [i for i in path_dict[sys][1]]
                     path_dict[newsys][1].append(sys)
 
-        visit(sys1)
+        visit(sys1.name)
 
         while True:
             min_distance = max(path_dict.values(), key=lambda x:x[0])[0]
@@ -478,10 +497,10 @@ def visit(sys):
                 break
             visit(newsys)
 
-        result = path_dict[sys2][1]
-        result.append(sys2)
+        result = path_dict[sys2.name][1]
+        result.append(sys2.name)
 
-        if result == [sys2]:
+        if result == [sys2.name]:
             raise KeyError("Two coordinate systems are not connected.")
         return result
 
diff --git a/sympy/diffgeom/rn.py b/sympy/diffgeom/rn.py
index 635e85f7365a..c6254734f4da 100644
--- a/sympy/diffgeom/rn.py
+++ b/sympy/diffgeom/rn.py
@@ -30,8 +30,8 @@
 r, theta = symbols('rho theta', nonnegative=True)
 
 relations_2d = {
-    ('rectangular', 'polar'): (sqrt(x**2 + y**2), atan2(y, x)),
-    ('polar', 'rectangular'): (r*cos(theta), r*sin(theta)),
+    ('rectangular', 'polar'): [(x, y), (sqrt(x**2 + y**2), atan2(y, x))],
+    ('polar', 'rectangular'): [(r, theta), (r*cos(theta), r*sin(theta))],
 }
 
 R2_r = CoordSystem('rectangular', R2_origin, (x, y), relations_2d)  # type: Any
@@ -74,18 +74,24 @@
 rho, psi, r, theta, phi = symbols('rho psi r theta phi', nonnegative=True)
 
 relations_3d = {
-    ('rectangular', 'cylindrical'): (sqrt(x**2 + y**2), atan2(y, x), z),
-    ('cylindrical', 'rectangular'): (rho*cos(psi), rho*sin(psi), z),
-    ('rectangular', 'spherical'): (sqrt(x**2 + y**2 + z**2),
-                                   acos(z/sqrt(x**2 + y**2 + z**2)),
-                                   atan2(y, x)),
-    ('spherical', 'rectangular'): (r*sin(theta)*cos(phi),
-                                   r*sin(theta)*sin(phi),
-                                   r*cos(theta)),
-    ('cylindrical', 'spherical'): (sqrt(rho**2 + z**2),
-                                   acos(z/sqrt(rho**2 + z**2)),
-                                   psi),
-    ('spherical', 'cylindrical'): (r*sin(theta), phi, r*cos(theta)),
+    ('rectangular', 'cylindrical'): [(x, y, z),
+                                     (sqrt(x**2 + y**2), atan2(y, x), z)],
+    ('cylindrical', 'rectangular'): [(rho, psi, z),
+                                     (rho*cos(psi), rho*sin(psi), z)],
+    ('rectangular', 'spherical'): [(x, y, z),
+                                   (sqrt(x**2 + y**2 + z**2),
+                                    acos(z/sqrt(x**2 + y**2 + z**2)),
+                                    atan2(y, x))],
+    ('spherical', 'rectangular'): [(r, theta, phi),
+                                   (r*sin(theta)*cos(phi),
+                                    r*sin(theta)*sin(phi),
+                                    r*cos(theta))],
+    ('cylindrical', 'spherical'): [(rho, psi, z),
+                                   (sqrt(rho**2 + z**2),
+                                    acos(z/sqrt(rho**2 + z**2)),
+                                    psi)],
+    ('spherical', 'cylindrical'): [(r, theta, phi),
+                                   (r*sin(theta), phi, r*cos(theta))],
 }
 
 R3_r = CoordSystem('rectangular', R3_origin, (x, y, z), relations_3d)  # type: Any
diff --git a/sympy/diffgeom/tests/test_diffgeom.py b/sympy/diffgeom/tests/test_diffgeom.py
index f5c1209702d4..dc36547e87eb 100644
--- a/sympy/diffgeom/tests/test_diffgeom.py
+++ b/sympy/diffgeom/tests/test_diffgeom.py
@@ -1,10 +1,10 @@
+from sympy.core import Lambda, Symbol, symbols
 from sympy.diffgeom.rn import R2, R2_p, R2_r, R3_r, R3_c, R3_s, R2_origin
 from sympy.diffgeom import (CoordSystem, Commutator, Differential, TensorProduct,
         WedgeProduct, BaseCovarDerivativeOp, CovarDerivativeOp, LieDerivative,
         covariant_order, contravariant_order, twoform_to_matrix, metric_to_Christoffel_1st,
         metric_to_Christoffel_2nd, metric_to_Riemann_components,
         metric_to_Ricci_components, intcurve_diffequ, intcurve_series)
-from sympy.core import Symbol, symbols
 from sympy.simplify import trigsimp, simplify
 from sympy.functions import sqrt, atan2, sin
 from sympy.matrices import Matrix
@@ -14,6 +14,78 @@
 TP = TensorProduct
 
 
+def test_coordsys_transform():
+    # test inverse transforms
+    p, q, r, s = symbols('p q r s')
+    rel = {('first', 'second'): [(p, q), (q, -p)]}
+    R2_pq = CoordSystem('first', R2_origin, [p, q], rel)
+    R2_rs = CoordSystem('second', R2_origin, [r, s], rel)
+    r, s = R2_rs.symbols
+    assert R2_rs.transform(R2_pq) == Matrix([[-s], [r]])
+
+    # inverse transform impossible case
+    a, b = symbols('a b', positive=True)
+    rel = {('first', 'second'): [(a,), (-a,)]}
+    R2_a = CoordSystem('first', R2_origin, [a], rel)
+    R2_b = CoordSystem('second', R2_origin, [b], rel)
+    # This transformation is uninvertible because there is no positive a, b satisfying a = -b
+    with raises(NotImplementedError):
+        R2_b.transform(R2_a)
+
+    # inverse transform ambiguous case
+    c, d = symbols('c d')
+    rel = {('first', 'second'): [(c,), (c**2,)]}
+    R2_c = CoordSystem('first', R2_origin, [c], rel)
+    R2_d = CoordSystem('second', R2_origin, [d], rel)
+    # The transform method should throw if it finds multiple inverses for a coordinate transformation.
+    with raises(ValueError):
+        R2_d.transform(R2_c)
+
+    # test indirect transformation
+    a, b, c, d, e, f = symbols('a, b, c, d, e, f')
+    rel = {('C1', 'C2'): [(a, b), (2*a, 3*b)],
+        ('C2', 'C3'): [(c, d), (3*c, 2*d)]}
+    C1 = CoordSystem('C1', R2_origin, (a, b), rel)
+    C2 = CoordSystem('C2', R2_origin, (c, d), rel)
+    C3 = CoordSystem('C3', R2_origin, (e, f), rel)
+    a, b = C1.symbols
+    c, d = C2.symbols
+    e, f = C3.symbols
+    assert C2.transform(C1) == Matrix([c/2, d/3])
+    assert C1.transform(C3) == Matrix([6*a, 6*b])
+    assert C3.transform(C1) == Matrix([e/6, f/6])
+    assert C3.transform(C2) == Matrix([e/3, f/2])
+
+    a, b, c, d, e, f = symbols('a, b, c, d, e, f')
+    rel = {('C1', 'C2'): [(a, b), (2*a, 3*b + 1)],
+        ('C3', 'C2'): [(e, f), (-e - 2, 2*f)]}
+    C1 = CoordSystem('C1', R2_origin, (a, b), rel)
+    C2 = CoordSystem('C2', R2_origin, (c, d), rel)
+    C3 = CoordSystem('C3', R2_origin, (e, f), rel)
+    a, b = C1.symbols
+    c, d = C2.symbols
+    e, f = C3.symbols
+    assert C2.transform(C1) == Matrix([c/2, (d - 1)/3])
+    assert C1.transform(C3) == Matrix([-2*a - 2, (3*b + 1)/2])
+    assert C3.transform(C1) == Matrix([-e/2 - 1, (2*f - 1)/3])
+    assert C3.transform(C2) == Matrix([-e - 2, 2*f])
+
+    # old signature uses Lambda
+    a, b, c, d, e, f = symbols('a, b, c, d, e, f')
+    rel = {('C1', 'C2'): Lambda((a, b), (2*a, 3*b + 1)),
+        ('C3', 'C2'): Lambda((e, f), (-e - 2, 2*f))}
+    C1 = CoordSystem('C1', R2_origin, (a, b), rel)
+    C2 = CoordSystem('C2', R2_origin, (c, d), rel)
+    C3 = CoordSystem('C3', R2_origin, (e, f), rel)
+    a, b = C1.symbols
+    c, d = C2.symbols
+    e, f = C3.symbols
+    assert C2.transform(C1) == Matrix([c/2, (d - 1)/3])
+    assert C1.transform(C3) == Matrix([-2*a - 2, (3*b + 1)/2])
+    assert C3.transform(C1) == Matrix([-e/2 - 1, (2*f - 1)/3])
+    assert C3.transform(C2) == Matrix([-e - 2, 2*f])
+
+
 def test_R2():
     x0, y0, r0, theta0 = symbols('x0, y0, r0, theta0', real=True)
     point_r = R2_r.point([x0, y0])
@@ -37,45 +109,6 @@ def test_R2():
             R2_r, R2_r.coord_tuple_transform_to(R2_p, m)).applyfunc(simplify)
 
 
-def test_inverse_transformations():
-    p, q, r, s = symbols('p q r s')
-
-    relations_quarter_rotation = {
-        ('first', 'second'): (q, -p)
-    }
-
-    R2_pq = CoordSystem('first', R2_origin, [p, q], relations_quarter_rotation)
-    R2_rs = CoordSystem('second', R2_origin, [r, s], relations_quarter_rotation)
-
-    # The transform method should derive the inverse transformation if not given explicitly
-    assert R2_rs.transform(R2_pq) == Matrix([[-R2_rs.symbols[1]], [R2_rs.symbols[0]]])
-
-    a, b = symbols('a b', positive=True)
-    relations_uninvertible_transformation = {
-        ('first', 'second'): (-a,)
-    }
-
-    R2_a = CoordSystem('first', R2_origin, [a], relations_uninvertible_transformation)
-    R2_b = CoordSystem('second', R2_origin, [b], relations_uninvertible_transformation)
-
-    # The transform method should throw if it cannot invert the coordinate transformation.
-    # This transformation is uninvertible because there is no positive a, b satisfying a = -b
-    with raises(NotImplementedError):
-        R2_b.transform(R2_a)
-
-    c, d = symbols('c d')
-    relations_ambiguous_inverse = {
-        ('first', 'second'): (c**2,)
-    }
-
-    R2_c = CoordSystem('first', R2_origin, [c], relations_ambiguous_inverse)
-    R2_d = CoordSystem('second', R2_origin, [d], relations_ambiguous_inverse)
-
-    # The transform method should throw if it finds multiple inverses for a coordinate transformation.
-    with raises(ValueError):
-        R2_d.transform(R2_c)
-
-
 def test_R3():
     a, b, c = symbols('a b c', positive=True)
     m = Matrix([[a], [b], [c]])


Code review comments:
diff_hunk: @@ -396,52 +401,66 @@ def transformation(self, sys):
         if self == sys:
             expr = Matrix(self.symbols)
         elif key in self.relations:
-            expr = Matrix(self.relations[key])
+            expr = Matrix(self.relations[key][1])
         elif key[::-1] in self.relations:
             expr = Matrix(self._inverse_transformation(sys, self))
         else:
             expr = Matrix(self._indirect_transformation(self, sys))
         return Lambda(signature, expr)
 
     @staticmethod
-    def _inverse_transformation(sys1, sys2):
-        # Find the transformation relation from sys2 to sys1
-        forward_transform_expressions = sys1.transform(sys2)
-
-        inv_results = solve(
-            [t[0] - t[1] for t in zip(sys2.symbols, forward_transform_expressions)],
-            list(sys1.symbols), dict=True)
-        if len(inv_results) == 0:
-            raise NotImplementedError(
-                "Cannot solve inverse of transformation from {} to {}".format(sys1, sys2))
-        elif len(inv_results) > 1:
-            raise ValueError(
-                "Obtained multiple results for inverse of transformation from {} to {}".format(sys1, sys2)
-            )
+    def _solve_inverse(sym1, sym2, exprs, sys1_name, sys2_name):
+        ret = solve(
+            [t[0] - t[1] for t in zip(sym2, exprs)],
+            list(sym1), dict=True)
+
+        if len(ret) == 0:
+            temp = "Cannot solve inverse relation from {} to {}."
+            raise NotImplementedError(temp.format(sys1_name, sys2_name))
+        elif len(ret) > 1:
+            temp = "Obtained multiple inverse relation from {} to {}."
+            raise ValueError(temp.format(sys1_name, sys2_name))
+
+        return ret[0]
 
-        inv_results = inv_results[0]
+    @classmethod
+    def _inverse_transformation(cls, sys1, sys2):
+        # Find the transformation relation from sys2 to sys1
+        forward = sys1.transform(sys2)
+        inv_results = cls._solve_inverse(sys1.symbols, sys2.symbols, forward,
+                                         sys1.name, sys2.name)
         signature = tuple(sys1.symbols)
         return [inv_results[s] for s in signature]
 
     @classmethod
     @cacheit
     def _indirect_transformation(cls, sys1, sys2):
-        # Find the transformation relation between two indirectly connected coordinate systems
+        # Find the transformation relation betweentwo indirectly connected
path: sympy/diffgeom/diffgeom.py
text: @user1:
```suggestion
        # Find the transformation relation between two indirectly connected
```