Repository: sympy/sympy
Pull request: define dummies once for geometry module
<!-- 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. -->


#### Brief description of what is fixed or changed

fixes #24581

#### 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 -->
* geometry
  * fixes memory leak caused by unnecessarily creating new Dummies
<!-- END RELEASE NOTES -->


Problem statement:
memory leak, not a cache issue
Hello

I have the following code that uses the latest sympy (1.11.1) tested with python 3.8.10 on Windows 10 and with python 3.11.1 on macOS 13.1, the memory is constantly growing. Any ideas why?

```
from sympy import Point, Ray, Circle
while True:
    circle = Circle(Point(0.0, 0.0), 0.5)
    ray = Ray(Point(0.2, 0.3), Point(0.3, 0.2))
    res = circle.intersection(ray)
```

NOTE: the growth in geometry has been addressed, but the scope of the issue is greater than that module so the issue is still open.

Hints from issue discussion:
I can reproduce this but no idea why it happens. Clearing the cache doesn't help. Calling `gc.collect` doesn't help.
With this diff the memory usage stays constant:
```diff
diff --git a/sympy/geometry/ellipse.py b/sympy/geometry/ellipse.py
index 0c1c5d0..35d92af 100644
--- a/sympy/geometry/ellipse.py
+++ b/sympy/geometry/ellipse.py
@@ -14,7 +14,7 @@
 from sympy.core.logic import fuzzy_bool
 from sympy.core.numbers import Rational, oo
 from sympy.core.sorting import ordered
-from sympy.core.symbol import Dummy, uniquely_named_symbol, _symbol
+from sympy.core.symbol import Symbol, Dummy, uniquely_named_symbol, _symbol
 from sympy.simplify import simplify, trigsimp
 from sympy.functions.elementary.miscellaneous import sqrt, Max
 from sympy.functions.elementary.trigonometric import cos, sin
@@ -664,8 +664,8 @@ def intersection(self, o):
         [Point2D(-17/5, -12/5), Point2D(-17/5, 12/5), Point2D(7/5, -12/5), Point2D(7/5, 12/5)]
         """
         # TODO: Replace solve with nonlinsolve, when nonlinsolve will be able to solve in real domain
-        x = Dummy('x', real=True)
-        y = Dummy('y', real=True)
+        x = Symbol('x', real=True)
+        y = Symbol('y', real=True)
 
         if isinstance(o, Point):
             if o in self:
```
So what is happening is that every time `intersection` is called two new `Dummy` symbols are created. Somewhere there is some sort of cache that grows without bound if more and more symbols are created.
It's here (and a similar one in fields.py):
https://github.com/sympy/sympy/blob/f8e33851e174bd686ac8dc88d75f955cf5dc8eeb/sympy/polys/rings.py#L194
This cache grows without bound every time a polynomial ring with new symbols is used.

A fix could be to use `lru_cache` instead of a global cache.

I am not entirely sure that this is only a cache though because it also associates with a particular PolyElement class and it is possible that that class needs to be globally unique.
Can't that be solved by using a WeakDictionary? 

I've never been a huge fan of this class design in the polys FWIW.
I think that a WeakValue dictionary could work but it would probably not achieve the performance benefits intended by the cache in the first place. The point is that this second call is faster:
```python
In [1]: %time K = QQ[x,y]
CPU times: user 4 ms, sys: 0 ns, total: 4 ms
Wall time: 2.38 ms

In [2]: %time K = QQ[x,y]
CPU times: user 0 ns, sys: 0 ns, total: 0 ns
Wall time: 89.6 µs
```
Usually these rings are created transiently as part of an operation like `factor`, `cancel` etc. Since the ring is only a temporary object there will probably not be many in existence at any one time so with a WeakValue dictionary the ring will always disappear. In that case we would be better off just not having any cache. On the other hand of many operations like `factor`, `cancel` etc are used as part of a complex operation it is likely that the same ring will be reconstructed many times so keeping a cache can save time on that.

I think that probably a LRU cache is good here but given the size of these rings it should probably be kept quite small:
```diff
diff --git a/sympy/polys/fields.py b/sympy/polys/fields.py
index a3f239c..cdb079a 100644
--- a/sympy/polys/fields.py
+++ b/sympy/polys/fields.py
@@ -6,6 +6,7 @@
 
 from operator import add, mul, lt, le, gt, ge
 
+from sympy.core.cache import __cacheit as _cacheit
 from sympy.core.expr import Expr
 from sympy.core.mod import Mod
 from sympy.core.numbers import Exp1
@@ -99,11 +100,11 @@ def sfield(exprs, *symbols, **options):
     else:
         return (_field, fracs)
 
-_field_cache: dict[Any, Any] = {}
 
 class FracField(DefaultPrinting):
     """Multivariate distributed rational function field. """
 
+    @_cacheit(10)
     def __new__(cls, symbols, domain, order=lex):
         from sympy.polys.rings import PolyRing
         ring = PolyRing(symbols, domain, order)
@@ -113,7 +114,7 @@ def __new__(cls, symbols, domain, order=lex):
         order = ring.order
 
         _hash_tuple = (cls.__name__, symbols, ngens, domain, order)
-        obj = _field_cache.get(_hash_tuple)
+        obj = None
 
         if obj is None:
             obj = object.__new__(cls)
@@ -138,8 +139,6 @@ def __new__(cls, symbols, domain, order=lex):
                     if not hasattr(obj, name):
                         setattr(obj, name, generator)
 
-            _field_cache[_hash_tuple] = obj
-
         return obj
 
     def _gens(self):
diff --git a/sympy/polys/rings.py b/sympy/polys/rings.py
index 0db1897..3a56486 100644
--- a/sympy/polys/rings.py
+++ b/sympy/polys/rings.py
@@ -7,6 +7,7 @@
 from functools import reduce
 from types import GeneratorType
 
+from sympy.core.cache import __cacheit as _cacheit
 from sympy.core.expr import Expr
 from sympy.core.numbers import igcd, oo
 from sympy.core.symbol import Symbol, symbols as _symbols
@@ -191,11 +192,11 @@ def _parse_symbols(symbols):
 
     raise GeneratorsError("expected a string, Symbol or expression or a non-empty sequence of strings, Symbols or expressions")
 
-_ring_cache: dict[Any, Any] = {}
 
 class PolyRing(DefaultPrinting, IPolys):
     """Multivariate distributed polynomial ring. """
 
+    @_cacheit(10)
     def __new__(cls, symbols, domain, order=lex):
         symbols = tuple(_parse_symbols(symbols))
         ngens = len(symbols)
@@ -203,7 +204,7 @@ def __new__(cls, symbols, domain, order=lex):
         order = OrderOpt.preprocess(order)
 
         _hash_tuple = (cls.__name__, symbols, ngens, domain, order)
-        obj = _ring_cache.get(_hash_tuple)
+        obj = None
 
         if obj is None:
             if domain.is_Composite and set(symbols) & set(domain.symbols):
@@ -257,8 +258,6 @@ def __new__(cls, symbols, domain, order=lex):
                     if not hasattr(obj, name):
                         setattr(obj, name, generator)
 
-            _ring_cache[_hash_tuple] = obj
-
         return obj
 
     def _gens(self):
```
In the long run the best solution would be to make the ring objects lighter weight and faster to construct.
> With this diff the memory usage stays constant:

Those two Dummy symbols should be created at the top of ellipse.py and then used the 4 times when needed in the various classes -- there is no need to keep creating new Dummy symbols.
To be clear plenty of other operations can cause this cache to grow and even to grow more quickly:
```
while True:
    cancel(Dummy('x') + 1)
```
I thought the point was that these objects break if they aren't singletonized? I seem to remember that sort of thing being the case for these classes in the polys, but correct me if I am wrong. If that's the case, you'd need a version of lru_cache that also acts like a weak dictionary (i.e., least recently used items aren't removed from the cache if they are also referenced somewhere). 
> Those two Dummy symbols should be created at the top of ellipse.py and then used the 4 times when needed in the various classes -- there is no need to keep creating new Dummy symbols.

Wouldn't that change the semantics? A dummy symbol is supposed to be unequal to everything except for itself. Two separate intersections should have dummies that are unequal, but reusing the same Dummy would make them equal. 
> I thought the point was that these objects break if they aren't singletonized?

Yes, actually they do. I just tried removing the cache and a bunch of poly tests failed. The problem is e.g. this:
https://github.com/sympy/sympy/blob/f8e33851e174bd686ac8dc88d75f955cf5dc8eeb/sympy/polys/rings.py#L419
Each ring dynamically creates a class for its elements and distinct copies of the same ring would have different classes.

I've opened #24585 which for now removes the ring cache altogether so we can see if there's a noticeable impact on performance and if the need for the ring cache is gone.

Longer term I would like to remove the dynamic class generation altogether. Also we should just make the construction of rings faster.

Commit message: define dummies once for module

Patch to review:
diff --git a/sympy/geometry/ellipse.py b/sympy/geometry/ellipse.py
index 0c1c5d06c4cb..bc38a5427e2f 100644
--- a/sympy/geometry/ellipse.py
+++ b/sympy/geometry/ellipse.py
@@ -34,6 +34,8 @@
 
 import random
 
+x = Dummy('x', real=True)
+y = Dummy('y', real=True)
 
 class Ellipse(GeometrySet):
     """An elliptical GeometryEntity.
@@ -103,9 +105,6 @@ class Ellipse(GeometrySet):
 
     def __contains__(self, o):
         if isinstance(o, Point):
-            x = Dummy('x', real=True)
-            y = Dummy('y', real=True)
-
             res = self.equation(x, y).subs({x: o.x, y: o.y})
             return trigsimp(simplify(res)) is S.Zero
         elif isinstance(o, Ellipse):
@@ -664,8 +663,6 @@ def intersection(self, o):
         [Point2D(-17/5, -12/5), Point2D(-17/5, 12/5), Point2D(7/5, -12/5), Point2D(7/5, 12/5)]
         """
         # TODO: Replace solve with nonlinsolve, when nonlinsolve will be able to solve in real domain
-        x = Dummy('x', real=True)
-        y = Dummy('y', real=True)
 
         if isinstance(o, Point):
             if o in self:
@@ -925,7 +922,6 @@ def normal_lines(self, p, prec=None):
 
         # find the 4 normal points and construct lines through them with
         # the corresponding slope
-        x, y = Dummy('x', real=True), Dummy('y', real=True)
         eq = self.equation(x, y)
         dydx = idiff(eq, y, x)
         norm = -1/dydx
@@ -1299,7 +1295,6 @@ def tangent_lines(self, p):
             # else p is outside the ellipse or we can't tell. In case of the
             # latter, the solutions returned will only be valid if
             # the point is not inside the ellipse; if it is, nan will result.
-            x, y = Dummy('x'), Dummy('y')
             eq = self.equation(x, y)
             dydx = idiff(eq, y, x)
             slope = Line(p, Point(x, y)).slope
diff --git a/sympy/geometry/entity.py b/sympy/geometry/entity.py
index 0a016bd22c40..f3d002537d40 100644
--- a/sympy/geometry/entity.py
+++ b/sympy/geometry/entity.py
@@ -64,6 +64,9 @@
 ]
 
 
+x, y = Dummy(), Dummy()
+T = Dummy('t', real=True)
+
 class GeometryEntity(Basic, EvalfMixin):
     """The base class for all geometrical entities.
 
@@ -392,15 +395,15 @@ def reflect(self, line):
         l = line
         o = Point(0, 0)
         if l.slope.is_zero:
-            y = l.args[0].y
-            if not y:  # x-axis
+            v = l.args[0].y
+            if not v:  # x-axis
                 return g.scale(y=-1)
-            reps = [(p, p.translate(y=2*(y - p.y))) for p in g.atoms(Point)]
+            reps = [(p, p.translate(y=2*(v - p.y))) for p in g.atoms(Point)]
         elif l.slope is oo:
-            x = l.args[0].x
-            if not x:  # y-axis
+            v = l.args[0].x
+            if not v:  # y-axis
                 return g.scale(x=-1)
-            reps = [(p, p.translate(x=2*(x - p.x))) for p in g.atoms(Point)]
+            reps = [(p, p.translate(x=2*(v - p.x))) for p in g.atoms(Point)]
         else:
             if not hasattr(g, 'reflect') and not all(
                     isinstance(arg, Point) for arg in g.args):
@@ -410,7 +413,6 @@ def reflect(self, line):
             c = l.coefficients
             d = -c[-1]/c[1]  # y-intercept
             # apply the transform to a single point
-            x, y = Dummy(), Dummy()
             xf = Point(x, y)
             xf = xf.translate(y=-d).rotate(-a, o).scale(y=-1
                 ).rotate(a, o).translate(y=d)
@@ -528,7 +530,6 @@ def parameter_value(self, other, t):
             other = Point(other, dim=self.ambient_dimension)
         if not isinstance(other, Point):
             raise ValueError("other must be a point")
-        T = Dummy('t', real=True)
         sol = solve(self.arbitrary_point(T) - other, T, dict=True)
         if not sol:
             raise ValueError("Given point is not on %s" % func_name(self))
diff --git a/sympy/geometry/line.py b/sympy/geometry/line.py
index b7c023ab1ad4..5a6911c3e272 100644
--- a/sympy/geometry/line.py
+++ b/sympy/geometry/line.py
@@ -45,6 +45,9 @@
 import random
 
 
+t, u = [Dummy(i) for i in 'tu']
+
+
 class LinearEntity(GeometrySet):
     """A base class for all linear entities (Line, Ray and Segment)
     in n-dimensional Euclidean space.
@@ -547,7 +550,6 @@ def intersect_parallel_segments(seg1, seg2):
                 # arbitrary points, when  equal, both give a
                 # non-negative parameter when the arbitrary point
                 # coordinates are equated
-                t, u = [Dummy(i) for i in 'tu']
                 tu = solve(self.arbitrary_point(t) - other.arbitrary_point(u),
                     t, u, dict=True)[0]
                 def ok(p, l):
@@ -1038,7 +1040,6 @@ def random_point(self, seed=None):
             rng = random.Random(seed)
         else:
             rng = random
-        t = Dummy()
         pt = self.arbitrary_point(t)
         if isinstance(self, Ray):
             v = abs(rng.gauss(0, 1))
diff --git a/sympy/geometry/plane.py b/sympy/geometry/plane.py
index 8c677bdc337d..f6bbdc8492ad 100644
--- a/sympy/geometry/plane.py
+++ b/sympy/geometry/plane.py
@@ -24,6 +24,9 @@
 import random
 
 
+x, y, z = map(Dummy, 'xyz')
+t = Dummy()  # intentionally name left blank
+
 class Plane(GeometryEntity):
     """
     A plane is a flat, two-dimensional surface. A plane is the two-dimensional
@@ -74,10 +77,8 @@ def __new__(cls, p1, a=None, b=None, **kwargs):
         return GeometryEntity.__new__(cls, p1, normal_vector, **kwargs)
 
     def __contains__(self, o):
-        x, y, z = map(Dummy, 'xyz')
         k = self.equation(x, y, z)
         if isinstance(o, (LinearEntity, LinearEntity3D)):
-            t = Dummy()
             d = Point3D(o.arbitrary_point(t))
             e = k.subs([(x, d.x), (y, d.y), (z, d.z)])
             return e.equals(0)
@@ -404,7 +405,6 @@ def intersection(self, o):
             if o in self:
                 return [o]
             else:
-                t = Dummy()  # unnamed else it may clash with a symbol in o
                 a = Point3D(o.arbitrary_point(t))
                 p1, n = self.p1, Point3D(self.normal_vector)
 
@@ -454,7 +454,6 @@ def is_coplanar(self, o):
         True
         """
         if isinstance(o, Plane):
-            x, y, z = map(Dummy, 'xyz')
             return not cancel(self.equation(x, y, z)/o.equation(x, y, z)).has(x, y, z)
         if isinstance(o, Point3D):
             return o in self
@@ -818,11 +817,10 @@ def random_point(self, seed=None):
             rng = random.Random(seed)
         else:
             rng = random
-        u, v = Dummy('u'), Dummy('v')
         params = {
-            u: 2*Rational(rng.gauss(0, 1)) - 1,
-            v: 2*Rational(rng.gauss(0, 1)) - 1}
-        return self.arbitrary_point(u, v).subs(params)
+            x: 2*Rational(rng.gauss(0, 1)) - 1,
+            y: 2*Rational(rng.gauss(0, 1)) - 1}
+        return self.arbitrary_point(x, y).subs(params)
 
     def parameter_value(self, other, u, v=None):
         """Return the parameter(s) corresponding to the given point.
diff --git a/sympy/geometry/polygon.py b/sympy/geometry/polygon.py
index 003568e9689b..4badbafb1f4b 100644
--- a/sympy/geometry/polygon.py
+++ b/sympy/geometry/polygon.py
@@ -22,6 +22,9 @@
 import warnings
 
 
+x, y, T = symbols('x, y, T', real=True, cls=Dummy)
+
+
 class Polygon(GeometrySet):
     """A two-dimensional polygon.
 
@@ -848,7 +851,6 @@ def parameter_value(self, other, t):
         if other.free_symbols:
             raise NotImplementedError('non-numeric coordinates')
         unknown = False
-        T = Dummy('t', real=True)
         p = self.arbitrary_point(T)
         for pt, cond in p.args:
             sol = solve(pt - other, T, dict=True)
@@ -1004,7 +1006,6 @@ def cut_section(self, line):
         points = list(self.vertices)
         points.append(points[0])
 
-        x, y = symbols('x, y', real=True, cls=Dummy)
         eq = line.equation(x, y)
 
         # considering equation of line to be `ax +by + c`


Code review comments:
diff_hunk: @@ -34,6 +34,8 @@
 
 import random
 
+x = Dummy('x', real=True)
+y = Dummy('y', real=True)
path: sympy/geometry/ellipse.py
text: @user1:
Would it be better to use more anonymous dummies `Dummy()` for anywhere where users seeing dummies should be a bug?
I think that it confuses, if such bug exist, at least in jupyter.

@author:
What about the last commit where a generic name is used. I have found debugging to be a little difficult at times when trying to track down a generic "dummy_123" or "x". Now if the symbol appears it will be obvious where it was generated.

@user2:
If a Dummy is returned to a user, it should use a mathematically meaningful name (typically a single letter name like x or t). I agree it's useful to name internal dummies for debugging. 

@author:
None of these should be returned.

@user2:
I guess I misunderstood. I thought these were dummies that are returned to the user as a new parameter in the solution. If they are only used internally, then I agree that they can be reused. For dummies like that, having a long (but also unique) name like `Dummy('geometry_internal')` would be best, so that it's very obvious where it came from if it ever does leak out.