Coverage for /usr/lib/python3/dist-packages/sympy/geometry/util.py: 9%

1"""Utility functions for geometrical entities.

3Contains

4========

5intersection

6convex_hull

7closest_points

8farthest_points

9are_coplanar

10are_similar

12"""

14from collections import deque

15from math import sqrt as _sqrt

18from .entity import GeometryEntity

19from .exceptions import GeometryError

20from .point import Point, Point2D, Point3D

21from sympy.core.containers import OrderedSet

22from sympy.core.exprtools import factor_terms

23from sympy.core.function import Function, expand_mul

24from sympy.core.sorting import ordered

25from sympy.core.symbol import Symbol

26from sympy.core.singleton import S

27from sympy.polys.polytools import cancel

28from sympy.functions.elementary.miscellaneous import sqrt

29from sympy.utilities.iterables import is_sequence

32def find(x, equation):

33 """

34 Checks whether a Symbol matching ``x`` is present in ``equation``

35 or not. If present, the matching symbol is returned, else a

36 ValueError is raised. If ``x`` is a string the matching symbol

37 will have the same name; if ``x`` is a Symbol then it will be

38 returned if found.

40 Examples

41 ========

43 >>> from sympy.geometry.util import find

44 >>> from sympy import Dummy

45 >>> from sympy.abc import x

46 >>> find('x', x)

47 x

48 >>> find('x', Dummy('x'))

49 _x

51 The dummy symbol is returned since it has a matching name:

53 >>> _.name == 'x'

54 True

55 >>> find(x, Dummy('x'))

56 Traceback (most recent call last):

57 ...

58 ValueError: could not find x

59 """

61 free = equation.free_symbols

62 xs = [i for i in free if (i.name if isinstance(x, str) else i) == x]

63 if not xs:

64 raise ValueError('could not find %s' % x)

65 if len(xs) != 1:

66 raise ValueError('ambiguous %s' % x)

67 return xs[0]

70def _ordered_points(p):

71 """Return the tuple of points sorted numerically according to args"""

72 return tuple(sorted(p, key=lambda x: x.args))

75def are_coplanar(*e):

76 """ Returns True if the given entities are coplanar otherwise False

78 Parameters

79 ==========

81 e: entities to be checked for being coplanar

83 Returns

84 =======

86 Boolean

88 Examples

89 ========

91 >>> from sympy import Point3D, Line3D

92 >>> from sympy.geometry.util import are_coplanar

93 >>> a = Line3D(Point3D(5, 0, 0), Point3D(1, -1, 1))

94 >>> b = Line3D(Point3D(0, -2, 0), Point3D(3, 1, 1))

95 >>> c = Line3D(Point3D(0, -1, 0), Point3D(5, -1, 9))

96 >>> are_coplanar(a, b, c)

97 False

99 """

100 from .line import LinearEntity3D

101 from .plane import Plane

102 # XXX update tests for coverage

103

104 e = set(e)

105 # first work with a Plane if present

106 for i in list(e):

107 if isinstance(i, Plane):

108 e.remove(i)

109 return all(p.is_coplanar(i) for p in e)

110

111 if all(isinstance(i, Point3D) for i in e):

112 if len(e) < 3:

113 return False

114

115 # remove pts that are collinear with 2 pts

116 a, b = e.pop(), e.pop()

117 for i in list(e):

118 if Point3D.are_collinear(a, b, i):

119 e.remove(i)

120

121 if not e:

122 return False

123 else:

124 # define a plane

125 p = Plane(a, b, e.pop())

126 for i in e:

127 if i not in p:

128 return False

129 return True

130 else:

131 pt3d = []

132 for i in e:

133 if isinstance(i, Point3D):

134 pt3d.append(i)

135 elif isinstance(i, LinearEntity3D):

136 pt3d.extend(i.args)

137 elif isinstance(i, GeometryEntity): # XXX we should have a GeometryEntity3D class so we can tell the difference between 2D and 3D -- here we just want to deal with 2D objects; if new 3D objects are encountered that we didn't handle above, an error should be raised

138 # all 2D objects have some Point that defines them; so convert those points to 3D pts by making z=0

139 for p in i.args:

140 if isinstance(p, Point):

141 pt3d.append(Point3D(*(p.args + (0,))))

142 return are_coplanar(*pt3d)

143

144

145def are_similar(e1, e2):

146 """Are two geometrical entities similar.

147

148 Can one geometrical entity be uniformly scaled to the other?

149

150 Parameters

151 ==========

152

153 e1 : GeometryEntity

154 e2 : GeometryEntity

155

156 Returns

157 =======

158

159 are_similar : boolean

160

161 Raises

162 ======

163

164 GeometryError

165 When `e1` and `e2` cannot be compared.

166

167 Notes

168 =====

169

170 If the two objects are equal then they are similar.

171

172 See Also

173 ========

174

175 sympy.geometry.entity.GeometryEntity.is_similar

176

177 Examples

178 ========

179

180 >>> from sympy import Point, Circle, Triangle, are_similar

181 >>> c1, c2 = Circle(Point(0, 0), 4), Circle(Point(1, 4), 3)

182 >>> t1 = Triangle(Point(0, 0), Point(1, 0), Point(0, 1))

183 >>> t2 = Triangle(Point(0, 0), Point(2, 0), Point(0, 2))

184 >>> t3 = Triangle(Point(0, 0), Point(3, 0), Point(0, 1))

185 >>> are_similar(t1, t2)

186 True

187 >>> are_similar(t1, t3)

188 False

189

190 """

191 if e1 == e2:

192 return True

193 is_similar1 = getattr(e1, 'is_similar', None)

194 if is_similar1:

195 return is_similar1(e2)

196 is_similar2 = getattr(e2, 'is_similar', None)

197 if is_similar2:

198 return is_similar2(e1)

199 n1 = e1.__class__.__name__

200 n2 = e2.__class__.__name__

201 raise GeometryError(

202 "Cannot test similarity between %s and %s" % (n1, n2))

203

204

205def centroid(*args):

206 """Find the centroid (center of mass) of the collection containing only Points,

207 Segments or Polygons. The centroid is the weighted average of the individual centroid

208 where the weights are the lengths (of segments) or areas (of polygons).

209 Overlapping regions will add to the weight of that region.

210

211 If there are no objects (or a mixture of objects) then None is returned.

212

213 See Also

214 ========

215

216 sympy.geometry.point.Point, sympy.geometry.line.Segment,

217 sympy.geometry.polygon.Polygon

218

219 Examples

220 ========

221

222 >>> from sympy import Point, Segment, Polygon

223 >>> from sympy.geometry.util import centroid

224 >>> p = Polygon((0, 0), (10, 0), (10, 10))

225 >>> q = p.translate(0, 20)

226 >>> p.centroid, q.centroid

227 (Point2D(20/3, 10/3), Point2D(20/3, 70/3))

228 >>> centroid(p, q)

229 Point2D(20/3, 40/3)

230 >>> p, q = Segment((0, 0), (2, 0)), Segment((0, 0), (2, 2))

231 >>> centroid(p, q)

232 Point2D(1, 2 - sqrt(2))

233 >>> centroid(Point(0, 0), Point(2, 0))

234 Point2D(1, 0)

235

236 Stacking 3 polygons on top of each other effectively triples the

237 weight of that polygon:

238

239 >>> p = Polygon((0, 0), (1, 0), (1, 1), (0, 1))

240 >>> q = Polygon((1, 0), (3, 0), (3, 1), (1, 1))

241 >>> centroid(p, q)

242 Point2D(3/2, 1/2)

243 >>> centroid(p, p, p, q) # centroid x-coord shifts left

244 Point2D(11/10, 1/2)

245

246 Stacking the squares vertically above and below p has the same

247 effect:

248

249 >>> centroid(p, p.translate(0, 1), p.translate(0, -1), q)

250 Point2D(11/10, 1/2)

251

252 """

253 from .line import Segment

254 from .polygon import Polygon

255 if args:

256 if all(isinstance(g, Point) for g in args):

257 c = Point(0, 0)

258 for g in args:

259 c += g

260 den = len(args)

261 elif all(isinstance(g, Segment) for g in args):

262 c = Point(0, 0)

263 L = 0

264 for g in args:

265 l = g.length

266 c += g.midpoint*l

267 L += l

268 den = L

269 elif all(isinstance(g, Polygon) for g in args):

270 c = Point(0, 0)

271 A = 0

272 for g in args:

273 a = g.area

274 c += g.centroid*a

275 A += a

276 den = A

277 c /= den

278 return c.func(*[i.simplify() for i in c.args])

279

280

281def closest_points(*args):

282 """Return the subset of points from a set of points that were

283 the closest to each other in the 2D plane.

284

285 Parameters

286 ==========

287

288 args

289 A collection of Points on 2D plane.

290

291 Notes

292 =====

293

294 This can only be performed on a set of points whose coordinates can

295 be ordered on the number line. If there are no ties then a single

296 pair of Points will be in the set.

297

298 Examples

299 ========

300

301 >>> from sympy import closest_points, Triangle

302 >>> Triangle(sss=(3, 4, 5)).args

303 (Point2D(0, 0), Point2D(3, 0), Point2D(3, 4))

304 >>> closest_points(*_)

305 {(Point2D(0, 0), Point2D(3, 0))}

306

307 References

308 ==========

309

310 .. [1] https://www.cs.mcgill.ca/~cs251/ClosestPair/ClosestPairPS.html

311

312 .. [2] Sweep line algorithm

313 https://en.wikipedia.org/wiki/Sweep_line_algorithm

314

315 """

316 p = [Point2D(i) for i in set(args)]

317 if len(p) < 2:

318 raise ValueError('At least 2 distinct points must be given.')

319

320 try:

321 p.sort(key=lambda x: x.args)

322 except TypeError:

323 raise ValueError("The points could not be sorted.")

324

325 if not all(i.is_Rational for j in p for i in j.args):

326 def hypot(x, y):

327 arg = x*x + y*y

328 if arg.is_Rational:

329 return _sqrt(arg)

330 return sqrt(arg)

331 else:

332 from math import hypot

333

334 rv = [(0, 1)]

335 best_dist = hypot(p[1].x - p[0].x, p[1].y - p[0].y)

336 i = 2

337 left = 0

338 box = deque([0, 1])

339 while i < len(p):

340 while left < i and p[i][0] - p[left][0] > best_dist:

341 box.popleft()

342 left += 1

343

344 for j in box:

345 d = hypot(p[i].x - p[j].x, p[i].y - p[j].y)

346 if d < best_dist:

347 rv = [(j, i)]

348 elif d == best_dist:

349 rv.append((j, i))

350 else:

351 continue

352 best_dist = d

353 box.append(i)

354 i += 1

355

356 return {tuple([p[i] for i in pair]) for pair in rv}

357

358

359def convex_hull(*args, polygon=True):

360 """The convex hull surrounding the Points contained in the list of entities.

361

362 Parameters

363 ==========

364

365 args : a collection of Points, Segments and/or Polygons

366

367 Optional parameters

368 ===================

369

370 polygon : Boolean. If True, returns a Polygon, if false a tuple, see below.

371 Default is True.

372

373 Returns

374 =======

375

376 convex_hull : Polygon if ``polygon`` is True else as a tuple `(U, L)` where

377 ``L`` and ``U`` are the lower and upper hulls, respectively.

378

379 Notes

380 =====

381

382 This can only be performed on a set of points whose coordinates can

383 be ordered on the number line.

384

385 See Also

386 ========

387

388 sympy.geometry.point.Point, sympy.geometry.polygon.Polygon

389

390 Examples

391 ========

392

393 >>> from sympy import convex_hull

394 >>> points = [(1, 1), (1, 2), (3, 1), (-5, 2), (15, 4)]

395 >>> convex_hull(*points)

396 Polygon(Point2D(-5, 2), Point2D(1, 1), Point2D(3, 1), Point2D(15, 4))

397 >>> convex_hull(*points, **dict(polygon=False))

398 ([Point2D(-5, 2), Point2D(15, 4)],

399 [Point2D(-5, 2), Point2D(1, 1), Point2D(3, 1), Point2D(15, 4)])

400

401 References

402 ==========

403

404 .. [1] https://en.wikipedia.org/wiki/Graham_scan

405

406 .. [2] Andrew's Monotone Chain Algorithm

407 (A.M. Andrew,

408 "Another Efficient Algorithm for Convex Hulls in Two Dimensions", 1979)

409 https://web.archive.org/web/20210511015444/http://geomalgorithms.com/a10-_hull-1.html

410

411 """

412 from .line import Segment

413 from .polygon import Polygon

414 p = OrderedSet()

415 for e in args:

416 if not isinstance(e, GeometryEntity):

417 try:

418 e = Point(e)

419 except NotImplementedError:

420 raise ValueError('%s is not a GeometryEntity and cannot be made into Point' % str(e))

421 if isinstance(e, Point):

422 p.add(e)

423 elif isinstance(e, Segment):

424 p.update(e.points)

425 elif isinstance(e, Polygon):

426 p.update(e.vertices)

427 else:

428 raise NotImplementedError(

429 'Convex hull for %s not implemented.' % type(e))

430

431 # make sure all our points are of the same dimension

432 if any(len(x) != 2 for x in p):

433 raise ValueError('Can only compute the convex hull in two dimensions')

434

435 p = list(p)

436 if len(p) == 1:

437 return p[0] if polygon else (p[0], None)

438 elif len(p) == 2:

439 s = Segment(p[0], p[1])

440 return s if polygon else (s, None)

441

442 def _orientation(p, q, r):

443 '''Return positive if p-q-r are clockwise, neg if ccw, zero if

444 collinear.'''

445 return (q.y - p.y)*(r.x - p.x) - (q.x - p.x)*(r.y - p.y)

446

447 # scan to find upper and lower convex hulls of a set of 2d points.

448 U = []

449 L = []

450 try:

451 p.sort(key=lambda x: x.args)

452 except TypeError:

453 raise ValueError("The points could not be sorted.")

454 for p_i in p:

455 while len(U) > 1 and _orientation(U[-2], U[-1], p_i) <= 0:

456 U.pop()

457 while len(L) > 1 and _orientation(L[-2], L[-1], p_i) >= 0:

458 L.pop()

459 U.append(p_i)

460 L.append(p_i)

461 U.reverse()

462 convexHull = tuple(L + U[1:-1])

463

464 if len(convexHull) == 2:

465 s = Segment(convexHull[0], convexHull[1])

466 return s if polygon else (s, None)

467 if polygon:

468 return Polygon(*convexHull)

469 else:

470 U.reverse()

471 return (U, L)

472

473def farthest_points(*args):

474 """Return the subset of points from a set of points that were

475 the furthest apart from each other in the 2D plane.

476

477 Parameters

478 ==========

479

480 args

481 A collection of Points on 2D plane.

482

483 Notes

484 =====

485

486 This can only be performed on a set of points whose coordinates can

487 be ordered on the number line. If there are no ties then a single

488 pair of Points will be in the set.

489

490 Examples

491 ========

492

493 >>> from sympy.geometry import farthest_points, Triangle

494 >>> Triangle(sss=(3, 4, 5)).args

495 (Point2D(0, 0), Point2D(3, 0), Point2D(3, 4))

496 >>> farthest_points(*_)

497 {(Point2D(0, 0), Point2D(3, 4))}

498

499 References

500 ==========

501

502 .. [1] https://code.activestate.com/recipes/117225-convex-hull-and-diameter-of-2d-point-sets/

503

504 .. [2] Rotating Callipers Technique

505 https://en.wikipedia.org/wiki/Rotating_calipers

506

507 """

508

509 def rotatingCalipers(Points):

510 U, L = convex_hull(*Points, **{"polygon": False})

511

512 if L is None:

513 if isinstance(U, Point):

514 raise ValueError('At least two distinct points must be given.')

515 yield U.args

516 else:

517 i = 0

518 j = len(L) - 1

519 while i < len(U) - 1 or j > 0:

520 yield U[i], L[j]

521 # if all the way through one side of hull, advance the other side

522 if i == len(U) - 1:

523 j -= 1

524 elif j == 0:

525 i += 1

526 # still points left on both lists, compare slopes of next hull edges

527 # being careful to avoid divide-by-zero in slope calculation

528 elif (U[i+1].y - U[i].y) * (L[j].x - L[j-1].x) > \

529 (L[j].y - L[j-1].y) * (U[i+1].x - U[i].x):

530 i += 1

531 else:

532 j -= 1

533

534 p = [Point2D(i) for i in set(args)]

535

536 if not all(i.is_Rational for j in p for i in j.args):

537 def hypot(x, y):

538 arg = x*x + y*y

539 if arg.is_Rational:

540 return _sqrt(arg)

541 return sqrt(arg)

542 else:

543 from math import hypot

544

545 rv = []

546 diam = 0

547 for pair in rotatingCalipers(args):

548 h, q = _ordered_points(pair)

549 d = hypot(h.x - q.x, h.y - q.y)

550 if d > diam:

551 rv = [(h, q)]

552 elif d == diam:

553 rv.append((h, q))

554 else:

555 continue

556 diam = d

557

558 return set(rv)

559

560

561def idiff(eq, y, x, n=1):

562 """Return ``dy/dx`` assuming that ``eq == 0``.

563

564 Parameters

565 ==========

566

567 y : the dependent variable or a list of dependent variables (with y first)

568 x : the variable that the derivative is being taken with respect to

569 n : the order of the derivative (default is 1)

570

571 Examples

572 ========

573

574 >>> from sympy.abc import x, y, a

575 >>> from sympy.geometry.util import idiff

576

577 >>> circ = x**2 + y**2 - 4

578 >>> idiff(circ, y, x)

579 -x/y

580 >>> idiff(circ, y, x, 2).simplify()

581 (-x**2 - y**2)/y**3

582

583 Here, ``a`` is assumed to be independent of ``x``:

584

585 >>> idiff(x + a + y, y, x)

586 -1

587

588 Now the x-dependence of ``a`` is made explicit by listing ``a`` after

589 ``y`` in a list.

590

591 >>> idiff(x + a + y, [y, a], x)

592 -Derivative(a, x) - 1

593

594 See Also

595 ========

596

597 sympy.core.function.Derivative: represents unevaluated derivatives

598 sympy.core.function.diff: explicitly differentiates wrt symbols

599

600 """

601 if is_sequence(y):

602 dep = set(y)

603 y = y[0]

604 elif isinstance(y, Symbol):

605 dep = {y}

606 elif isinstance(y, Function):

607 pass

608 else:

609 raise ValueError("expecting x-dependent symbol(s) or function(s) but got: %s" % y)

610

611 f = {s: Function(s.name)(x) for s in eq.free_symbols

612 if s != x and s in dep}

613

614 if isinstance(y, Symbol):

615 dydx = Function(y.name)(x).diff(x)

616 else:

617 dydx = y.diff(x)

618

619 eq = eq.subs(f)

620 derivs = {}

621 for i in range(n):

622 # equation will be linear in dydx, a*dydx + b, so dydx = -b/a

623 deq = eq.diff(x)

624 b = deq.xreplace({dydx: S.Zero})

625 a = (deq - b).xreplace({dydx: S.One})

626 yp = factor_terms(expand_mul(cancel((-b/a).subs(derivs)), deep=False))

627 if i == n - 1:

628 return yp.subs([(v, k) for k, v in f.items()])

629 derivs[dydx] = yp

630 eq = dydx - yp

631 dydx = dydx.diff(x)

632

633

634def intersection(*entities, pairwise=False, **kwargs):

635 """The intersection of a collection of GeometryEntity instances.

636

637 Parameters

638 ==========

639 entities : sequence of GeometryEntity

640 pairwise (keyword argument) : Can be either True or False

641

642 Returns

643 =======

644 intersection : list of GeometryEntity

645

646 Raises

647 ======

648 NotImplementedError

649 When unable to calculate intersection.

650

651 Notes

652 =====

653 The intersection of any geometrical entity with itself should return

654 a list with one item: the entity in question.

655 An intersection requires two or more entities. If only a single

656 entity is given then the function will return an empty list.

657 It is possible for `intersection` to miss intersections that one

658 knows exists because the required quantities were not fully

659 simplified internally.

660 Reals should be converted to Rationals, e.g. Rational(str(real_num))

661 or else failures due to floating point issues may result.

662

663 Case 1: When the keyword argument 'pairwise' is False (default value):

664 In this case, the function returns a list of intersections common to

665 all entities.

666

667 Case 2: When the keyword argument 'pairwise' is True:

668 In this case, the functions returns a list intersections that occur

669 between any pair of entities.

670

671 See Also

672 ========

673

674 sympy.geometry.entity.GeometryEntity.intersection

675

676 Examples

677 ========

678

679 >>> from sympy import Ray, Circle, intersection

680 >>> c = Circle((0, 1), 1)

681 >>> intersection(c, c.center)

682 []

683 >>> right = Ray((0, 0), (1, 0))

684 >>> up = Ray((0, 0), (0, 1))

685 >>> intersection(c, right, up)

686 [Point2D(0, 0)]

687 >>> intersection(c, right, up, pairwise=True)

688 [Point2D(0, 0), Point2D(0, 2)]

689 >>> left = Ray((1, 0), (0, 0))

690 >>> intersection(right, left)

691 [Segment2D(Point2D(0, 0), Point2D(1, 0))]

692

693 """

694 if len(entities) <= 1:

695 return []

696

697 # entities may be an immutable tuple

698 entities = list(entities)

699 for i, e in enumerate(entities):

700 if not isinstance(e, GeometryEntity):

701 entities[i] = Point(e)

702

703 if not pairwise:

704 # find the intersection common to all objects

705 res = entities[0].intersection(entities[1])

706 for entity in entities[2:]:

707 newres = []

708 for x in res:

709 newres.extend(x.intersection(entity))

710 res = newres

711 return res

712

713 # find all pairwise intersections

714 ans = []

715 for j in range(len(entities)):

716 for k in range(j + 1, len(entities)):

717 ans.extend(intersection(entities[j], entities[k]))

718 return list(ordered(set(ans)))