Repository: sympy/sympy

Issue:
implement parameter finding routine
#13539 started in an interesting direction by giving a parameter-finding function. This could be extended to solve for the parameter of any geometric object that implements the `arbitrary_point` method:

```diff
diff --git a/sympy/geometry/entity.py b/sympy/geometry/entity.py
index ef9f3d0..dd152f3 100644
--- a/sympy/geometry/entity.py
+++ b/sympy/geometry/entity.py
@@ -496,6 +496,37 @@ def translate(self, x=0, y=0):
                 newargs.append(a)
         return self.func(*newargs)

+    def parameter(self, other):
+        """Return the parameter corresponding to the given point. Evaluating
+        the entity at this parameter value will return the given point.
+
+        Examples
+        ========
+
+        >>> from sympy import Line, Point
+        >>> from sympy.abc import t
+        >>> a = Point(0,0)
+        >>> b = Point(2, 2)
+        >>> Line(a, b).parameter((1, 1))
+        1/2
+        >>> Line(a, b).arbitrary_point(t).subs(t, _)
+        Point2D(1, 1)
+        """
+        from sympy.solvers.solvers import solve
+        from sympy.geometry.point import Point
+        from sympy.core.symbol import Dummy
+        if not isinstance(other, GeometryEntity):
+            other = Point(other, dim=self.ambient_dimension)
+        if not isinstance(other, Point):
+            raise ValueError('other must be a Point')
+        k = Dummy()
+        sol = solve(self.arbitrary_point(k) - other, dict=True) or [{}]
+        v = sol[0].get(k, None)
+        if v is None:
+            raise ValueError('given point is not on the entity')
+        return v
+
+
 class GeometrySet(GeometryEntity, Set):
     """Parent class of all GeometryEntity that are also Sets
     (compatible with sympy.sets)
```

Tests for the different object would be added to the test suite (similar to the test given in the doctest above) for those objects that implement `arbitrary_point`.

Discussion:
smichr: @normalhuman , could you please take a look at this, especially the Plane methods `random_point` and `arbitrary_point`?
smichr: ping @saurabhkgp21 , since you worked on the first part, perhaps you could take a look at the proposed additions.

Patch:
diff --git a/sympy/geometry/curve.py b/sympy/geometry/curve.py
index 4cae33ba04..e8d5ecc37b 100644
--- a/sympy/geometry/curve.py
+++ b/sympy/geometry/curve.py
@@ -165,6 +165,10 @@ def free_symbols(self):
         free = free.difference({self.parameter})
         return free
 
+    @property
+    def ambient_dimension(self):
+        return len(self.args[0])
+
     @property
     def functions(self):
         """The functions specifying the curve.
diff --git a/sympy/geometry/entity.py b/sympy/geometry/entity.py
index d0bb07a7db..83c2e61822 100644
--- a/sympy/geometry/entity.py
+++ b/sympy/geometry/entity.py
@@ -25,10 +25,12 @@
 from sympy.core.compatibility import is_sequence
 from sympy.core.containers import Tuple
 from sympy.core.basic import Basic
+from sympy.core.symbol import _symbol
 from sympy.core.sympify import sympify
 from sympy.functions import cos, sin
 from sympy.matrices import eye
 from sympy.sets import Set
+from sympy.utilities.misc import func_name
 
 # How entities are ordered; used by __cmp__ in GeometryEntity
 ordering_of_classes = [
@@ -456,7 +458,7 @@ def scale(self, x=1, y=1, pt=None):
         Triangle(Point2D(1, 0), Point2D(-1/2, sqrt(3)/2), Point2D(-1/2, -sqrt(3)/2))
         >>> t.scale(2)
         Triangle(Point2D(2, 0), Point2D(-1, sqrt(3)/2), Point2D(-1, -sqrt(3)/2))
-        >>> t.scale(2,2)
+        >>> t.scale(2, 2)
         Triangle(Point2D(2, 0), Point2D(-1, sqrt(3)), Point2D(-1, -sqrt(3)))
 
         """
@@ -496,6 +498,37 @@ def translate(self, x=0, y=0):
                 newargs.append(a)
         return self.func(*newargs)
 
+    def parameter_value(self, other, t):
+        """Return the parameter corresponding to the given point.
+        Evaluating an arbitrary point of the entity at this parameter
+        value will return the given point.
+
+        Examples
+        ========
+
+        >>> from sympy import Line, Point
+        >>> from sympy.abc import t
+        >>> a = Point(0, 0)
+        >>> b = Point(2, 2)
+        >>> Line(a, b).parameter_value((1, 1), t)
+        {t: 1/2}
+        >>> Line(a, b).arbitrary_point(t).subs(_)
+        Point2D(1, 1)
+        """
+        from sympy.geometry.point import Point
+        from sympy.core.symbol import Dummy
+        from sympy.solvers.solvers import solve
+        if not isinstance(other, GeometryEntity):
+            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))
+        return {t: sol[0][T]}
+
+
 class GeometrySet(GeometryEntity, Set):
     """Parent class of all GeometryEntity that are also Sets
     (compatible with sympy.sets)
diff --git a/sympy/geometry/line.py b/sympy/geometry/line.py
index 485f721ad4..447cbdb8a2 100644
--- a/sympy/geometry/line.py
+++ b/sympy/geometry/line.py
@@ -22,8 +22,9 @@
 
 from sympy.core import S, sympify
 from sympy.core.compatibility import ordered
+from sympy.core.numbers import Rational
 from sympy.core.relational import Eq
-from sympy.core.symbol import _symbol
+from sympy.core.symbol import _symbol, Dummy
 from sympy.functions.elementary.trigonometric import (_pi_coeff as pi_coeff, acos, tan, atan2)
 from sympy.functions.elementary.piecewise import Piecewise
 from sympy.logic.boolalg import And
@@ -978,7 +979,7 @@ def proj_point(p):
         raise GeometryError(
             "Do not know how to project %s onto %s" % (other, self))
 
-    def random_point(self):
+    def random_point(self, seed=None):
         """A random point on a LinearEntity.
 
         Returns
@@ -994,32 +995,37 @@ def random_point(self):
         Examples
         ========
 
-        >>> from sympy import Point, Line
+        >>> from sympy import Point, Line, Ray, Segment
         >>> p1, p2 = Point(0, 0), Point(5, 3)
-        >>> l1 = Line(p1, p2)
-        >>> p3 = l1.random_point()
-        >>> # random point - don't know its coords in advance
-        >>> p3 # doctest: +ELLIPSIS
-        Point2D(...)
-        >>> # point should belong to the line
-        >>> p3 in l1
+        >>> line = Line(p1, p2)
+        >>> r = line.random_point(seed=42)  # seed value is optional
+        >>> r.n(3)
+        Point2D(-0.72, -0.432)
+        >>> r in line
         True
+        >>> Ray(p1, p2).random_point(seed=42).n(3)
+        Point2D(0.72, 0.432)
+        >>> Segment(p1, p2).random_point(seed=42).n(3)
+        Point2D(3.2, 1.92)
 
         """
-        from random import randint
-        from sympy.functions import floor
-
-        # The lower and upper
-        lower, upper = -2**32 - 1, 2**32
+        import random
 
+        if seed is not None:
+            rng = random.Random(seed)
+        else:
+            rng = random
+        t = Dummy()
+        pt = self.arbitrary_point(t)
         if isinstance(self, Ray):
-            lower = 0
-        if isinstance(self, Segment):
-            lower = 0
-            upper = floor(self.length)
-        t = randint(lower, upper)
-
-        return self.direction*t/abs(self.direction) + self.p1
+            v = abs(rng.gauss(0, 1))
+        elif isinstance(self, Segment):
+            v = rng.random()
+        elif isinstance(self, Line):
+            v = rng.gauss(0, 1)
+        else:
+            raise NotImplementedError('unhandled line type')
+        return pt.subs(t, Rational(v))
 
 
 class Line(LinearEntity):
diff --git a/sympy/geometry/plane.py b/sympy/geometry/plane.py
index e3c2be83c9..85c0f6375d 100644
--- a/sympy/geometry/plane.py
+++ b/sympy/geometry/plane.py
@@ -9,13 +9,14 @@
 
 from sympy import simplify
 from sympy.core import Dummy, Rational, S, Symbol
+from sympy.core.symbol import _symbol
 from sympy.core.compatibility import is_sequence
 from sympy.functions.elementary.trigonometric import cos, sin, acos, asin, sqrt
 from sympy.matrices import Matrix
 from sympy.polys.polytools import cancel
 from sympy.solvers import solve, linsolve
-from sympy.utilities.misc import filldedent
 from sympy.utilities.iterables import uniq
+from sympy.utilities.misc import filldedent, func_name
 
 from .entity import GeometryEntity
 from .point import Point, Point3D
@@ -136,21 +137,30 @@ def angle_between(self, o):
 
 
     def arbitrary_point(self, u=None, v=None):
-        """ Returns an arbitrary point on the Plane. If given two parameters,
-        the point ranges over the entire plane. If given one or no parameters,
-        returns a point with one parameter, varying which from 0 to 2*pi will
-        move the point in a circle of radius 1 about p1 of the Plane.
+        """ Returns an arbitrary point on the Plane. If given two
+        parameters, the point ranges over the entire plane. If given 1
+        or no parameters, returns a point with one parameter which,
+        when varying from 0 to 2*pi, moves the point in a circle of
+        radius 1 about p1 of the Plane.
 
         Examples
         ========
 
-        >>> from sympy.geometry.plane import Plane
-        >>> from sympy.abc import u, v, t
-        >>> p = Plane((3, 1, 4), (2, 0, 1), (-1, -2, -2))
+        >>> from sympy.geometry import Plane, Ray
+        >>> from sympy.abc import u, v, t, r
+        >>> p = Plane((1, 1, 1), normal_vector=(1, 0, 0))
         >>> p.arbitrary_point(u, v)
-        Point3D(-6*u - 3*v + 3, -3*u + 6*v + 1, 45*v + 4)
+        Point3D(1, u + 1, v + 1)
         >>> p.arbitrary_point(t)
-        Point3D(-sqrt(230)*sin(t)/230 - 2*sqrt(5)*cos(t)/5 + 3, sqrt(230)*sin(t)/115 - sqrt(5)*cos(t)/5 + 1, 3*sqrt(230)*sin(t)/46 + 4)
+        Point3D(1, cos(t) + 1, sin(t) + 1)
+
+        While arbitrary values of u and v can move the point anywhere in
+        the plane, the single-parameter point can be used to construct a
+        ray whose arbitrary point can be located at angle t and radius
+        r from p.p1:
+
+        >>> Ray(p.p1, _).arbitrary_point(r)
+        Point3D(1, r*cos(t) + 1, r*sin(t) + 1)
 
         Returns
         =======
@@ -159,8 +169,11 @@ def arbitrary_point(self, u=None, v=None):
 
         """
         circle = v is None
-        u = u or Dummy('u', real=True)
-        v = v or Dummy('v', real=True)
+        if circle:
+            u = _symbol(u or 't', real=True)
+        else:
+            u = _symbol(u or 'u', real=True)
+            v = _symbol(v or 'v', real=True)
         x, y, z = self.normal_vector
         a, b, c = self.p1.args
         # x1, y1, z1 is a nonzero vector parallel to the plane
@@ -777,7 +790,6 @@ def projection(self, pt):
             return rv
         return self.intersection(Line3D(rv, rv + Point3D(self.normal_vector)))[0]
 
-
     def random_point(self, seed=None):
         """ Returns a random point on the Plane.
 
@@ -786,6 +798,21 @@ def random_point(self, seed=None):
 
         Point3D
 
+        Examples
+        ========
+
+        >>> from sympy import Plane
+        >>> p = Plane((1, 0, 0), normal_vector=(0, 1, 0))
+        >>> r = p.random_point(seed=42)  # seed value is optional
+        >>> r.n(3)
+        Point3D(2.29, 0, -1.35)
+
+        The random point can be moved to lie on the circle of radius
+        1 centered on p1:
+
+        >>> c = p.p1 + (r - p.p1).unit
+        >>> c.distance(p.p1).equals(1)
+        True
         """
         import random
         if seed is not None:
@@ -793,5 +820,71 @@ def random_point(self, seed=None):
         else:
             rng = random
         u, v = Dummy('u'), Dummy('v')
-        params = {u: Rational(rng.random()), v: Rational(rng.random())}
+        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)
+
+    def parameter_value(self, other, u, v=None):
+        """Return the parameter(s) corresponding to the given point.
+
+        Examples
+        ========
+
+        >>> from sympy import Plane, Point, pi
+        >>> from sympy.abc import t, u, v
+        >>> p = Plane((2, 0, 0), (0, 0, 1), (0, 1, 0))
+
+        By default, the parameter value returned defines a point
+        that is a distance of 1 from the Plane's p1 value and
+        in line with the given point:
+
+        >>> on_circle = p.arbitrary_point(t).subs(t, pi/4)
+        >>> on_circle.distance(p.p1)
+        1
+        >>> p.parameter_value(on_circle, t)
+        {t: pi/4}
+
+        Moving the point twice as far from p1 does not change
+        the parameter value:
+
+        >>> off_circle = p.p1 + (on_circle - p.p1)*2
+        >>> off_circle.distance(p.p1)
+        2
+        >>> p.parameter_value(off_circle, t)
+        {t: pi/4}
+
+        If the 2-value parameter is desired, supply the two
+        parameter symbols and a replacement dictionary will
+        be returned:
+
+        >>> p.parameter_value(on_circle, u, v)
+        {u: sqrt(10)/10, v: sqrt(10)/30}
+        >>> p.parameter_value(off_circle, u, v)
+        {u: sqrt(10)/5, v: sqrt(10)/15}
+        """
+        from sympy.geometry.point import Point
+        from sympy.core.symbol import Dummy
+        from sympy.solvers.solvers import solve
+        if not isinstance(other, GeometryEntity):
+            other = Point(other, dim=self.ambient_dimension)
+        if not isinstance(other, Point):
+            raise ValueError("other must be a point")
+        if other == self.p1:
+            return other
+        if isinstance(u, Symbol) and v is None:
+            delta = self.arbitrary_point(u) - self.p1
+            eq = delta - (other - self.p1).unit
+            sol = solve(eq, u, dict=True)
+        elif isinstance(u, Symbol) and isinstance(v, Symbol):
+            pt = self.arbitrary_point(u, v)
+            sol = solve(pt - other, (u, v), dict=True)
+        else:
+            raise ValueError('expecting 1 or 2 symbols')
+        if not sol:
+            raise ValueError("Given point is not on %s" % func_name(self))
+        return sol[0]  # {t: tval} or {u: uval, v: vval}
+
+    @property
+    def ambient_dimension(self):
+        return self.p1.ambient_dimension
diff --git a/sympy/geometry/polygon.py b/sympy/geometry/polygon.py
index a764280f4c..d9cd2c7005 100644
--- a/sympy/geometry/polygon.py
+++ b/sympy/geometry/polygon.py
@@ -2,7 +2,7 @@
 
 from sympy.core import Expr, S, Symbol, oo, pi, sympify
 from sympy.core.compatibility import as_int, range, ordered
-from sympy.core.symbol import _symbol
+from sympy.core.symbol import _symbol, Dummy
 from sympy.functions.elementary.complexes import sign
 from sympy.functions.elementary.piecewise import Piecewise
 from sympy.functions.elementary.trigonometric import cos, sin, tan
@@ -12,6 +12,7 @@
 from sympy.simplify import simplify
 from sympy.utilities import default_sort_key
 from sympy.utilities.iterables import has_dups, has_variety, uniq
+from sympy.utilities.misc import func_name
 
 from .entity import GeometryEntity, GeometrySet
 from .point import Point
@@ -269,6 +270,10 @@ def angles(self):
                 ret[b] = ang
         return ret
 
+    @property
+    def ambient_dimension(self):
+        return self.vertices[0].ambient_dimension
+
     @property
     def perimeter(self):
         """The perimeter of the polygon.
@@ -607,6 +612,29 @@ def arbitrary_point(self, parameter='t'):
             perim_fraction_start = perim_fraction_end
         return Piecewise(*sides)
 
+    def parameter_value(self, other, t):
+        from sympy.solvers.solvers import solve
+        if not isinstance(other,GeometryEntity):
+            other = Point(other, dim=self.ambient_dimension)
+        if not isinstance(other,Point):
+            raise ValueError("other must be a point")
+        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)
+            if not sol:
+                continue
+            value = sol[0][T]
+            if simplify(cond.subs(T, value)) == True:
+                return {t: value}
+            unknown = True
+        if unknown:
+            raise ValueError("Given point may not be on %s" % func_name(self))
+        raise ValueError("Given point is not on %s" % func_name(self))
+
     def plot_interval(self, parameter='t'):
         """The plot interval for the default geometric plot of the polygon.
 


Test patch:
diff --git a/sympy/geometry/tests/test_curve.py b/sympy/geometry/tests/test_curve.py
index 1dd84312f1..fafb40e533 100644
--- a/sympy/geometry/tests/test_curve.py
+++ b/sympy/geometry/tests/test_curve.py
@@ -99,3 +99,10 @@ def test_length():
 
     c3 = Curve((t ** 2, t), (t, 2, 5))
     assert c3.length == -sqrt(17) - asinh(4) / 4 + asinh(10) / 4 + 5 * sqrt(101) / 2
+
+
+def test_parameter_value():
+    t = Symbol('t')
+    C = Curve([2*t, t**2], (t, 0, 2))
+    assert C.parameter_value((2, 1), t) == {t: 1}
+    raises(ValueError, lambda: C.parameter_value((2, 0), t))
diff --git a/sympy/geometry/tests/test_ellipse.py b/sympy/geometry/tests/test_ellipse.py
index 2dec7be2cd..4efd449c28 100644
--- a/sympy/geometry/tests/test_ellipse.py
+++ b/sympy/geometry/tests/test_ellipse.py
@@ -194,7 +194,7 @@ def line_list_close(ll1, ll2, prec):
     assert e1.intersection(Circle(Point(0, 2), 1)) == [Point(0, 1)]
     assert e1.intersection(Circle(Point(5, 0), 1)) == []
     assert e1.intersection(Ellipse(Point(2, 0), 1, 1)) == [Point(1, 0)]
-    assert e1.intersection(Ellipse(Point(5, 0), 1, 1,)) == []
+    assert e1.intersection(Ellipse(Point(5, 0), 1, 1)) == []
     assert e1.intersection(Point(2, 0)) == []
     assert e1.intersection(e1) == e1
     assert intersection(Ellipse(Point(0, 0), 2, 1), Ellipse(Point(3, 0), 1, 2)) == [Point(2, 0)]
@@ -205,7 +205,7 @@ def line_list_close(ll1, ll2, prec):
     assert Circle((0, 0), 1/2).intersection(
         Triangle((-1, 0), (1, 0), (0, 1))) == [
         Point(-1/2, 0), Point(1/2, 0)]
-    raises(TypeError, lambda: intersection(e2, Line((0, 0, 0), (0,0,1))))
+    raises(TypeError, lambda: intersection(e2, Line((0, 0, 0), (0, 0, 1))))
     raises(TypeError, lambda: intersection(e2, Rational(12)))
     # some special case intersections
     csmall = Circle(p1, 3)
@@ -325,7 +325,7 @@ def test_transform():
 
 
 def test_bounds():
-    e1 = Ellipse(Point(0,0), 3, 5)
+    e1 = Ellipse(Point(0, 0), 3, 5)
     e2 = Ellipse(Point(2, -2), 7, 7)
     c1 = Circle(Point(2, -2), 7)
     c2 = Circle(Point(-2, 0), Point(0, 2), Point(2, 0))
@@ -348,7 +348,7 @@ def test_reflect():
 
 
 def test_is_tangent():
-    e1 = Ellipse(Point(0,0), 3, 5)
+    e1 = Ellipse(Point(0, 0), 3, 5)
     c1 = Circle(Point(2, -2), 7)
     assert e1.is_tangent(Point(0, 0)) is False
     assert e1.is_tangent(Point(3, 0)) is False
@@ -385,3 +385,10 @@ def test_is_tangent():
     assert e1.is_tangent(Polygon((3, 0), (5, 7), (6, -5))) is False
     raises(TypeError, lambda: e1.is_tangent(Point(0, 0, 0)))
     raises(TypeError, lambda: e1.is_tangent(Rational(5)))
+
+
+def test_parameter_value():
+    t = Symbol('t')
+    e = Ellipse(Point(0, 0), 3, 5)
+    assert e.parameter_value((3, 0), t) == {t: 0}
+    raises(ValueError, lambda: e.parameter_value((4, 0), t))
diff --git a/sympy/geometry/tests/test_line.py b/sympy/geometry/tests/test_line.py
index 343216a64e..54b8cb350a 100644
--- a/sympy/geometry/tests/test_line.py
+++ b/sympy/geometry/tests/test_line.py
@@ -715,3 +715,11 @@ def _check():
     # midline intersection
     c, d = (-2, -3), (-2, 0)
     _check()
+
+
+def test_parameter_value():
+    t = Symbol('t')
+    p1, p2 = Point(0, 1), Point(5, 6)
+    l = Line(p1, p2)
+    assert l.parameter_value((5, 6), t) == {t: 1}
+    raises(ValueError, lambda: l.parameter_value((0, 0), t))
diff --git a/sympy/geometry/tests/test_plane.py b/sympy/geometry/tests/test_plane.py
index 37e297c7ab..11e9713846 100644
--- a/sympy/geometry/tests/test_plane.py
+++ b/sympy/geometry/tests/test_plane.py
@@ -227,3 +227,11 @@ def test_dimension_normalization():
     assert p.perpendicular_plane(a, b) == Plane(Point3D(0, 0, 0), (1, 0, 0))
     assert Plane((1, 2, 1), (2, 1, 0), (3, 1, 2)
         ).intersection((2, 1)) == [Point(2, 1, 0)]
+
+
+def test_parameter_value():
+    t, u, v = symbols("t, u v")
+    p = Plane((0, 0, 0), (0, 0, 1), (0, 1, 0))
+    assert p.parameter_value((0, -3, 2), t) == {t: asin(2*sqrt(13)/13)}
+    assert p.parameter_value((0, -3, 2), u, v) == {u: 3, v: 2}
+    raises(ValueError, lambda: p.parameter_value((1, 0, 0), t))
diff --git a/sympy/geometry/tests/test_polygon.py b/sympy/geometry/tests/test_polygon.py
index 3a60d2f0f5..f77ea78208 100644
--- a/sympy/geometry/tests/test_polygon.py
+++ b/sympy/geometry/tests/test_polygon.py
@@ -411,9 +411,18 @@ def test_intersection():
     assert poly1.intersection(RegularPolygon((-12, -15), 3, 3)) == []
 
 
+def test_parameter_value():
+    t = Symbol('t')
+    sq = Polygon((0, 0), (0, 1), (1, 1), (1, 0))
+    assert sq.parameter_value((0.5, 1), t) == {t: 3/8}
+    q = Polygon((0, 0), (2, 1), (2, 4), (4, 0))
+    assert q.parameter_value((4, 0), t) == {t: -6 + 3*sqrt(5)}  # ~= 0.708
+    raises(ValueError, lambda: sq.parameter_value((5, 6), t))
+
+
 def test_issue_12966():
     poly = Polygon(Point(0, 0), Point(0, 10), Point(5, 10), Point(5, 5),
-    Point(10, 5), Point(10, 0))
+        Point(10, 5), Point(10, 0))
     t = Symbol('t')
     pt = poly.arbitrary_point(t)
     DELTA = 5/poly.perimeter