plot_implicit.py 15 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431
  1. """Implicit plotting module for SymPy.
  2. Explanation
  3. ===========
  4. The module implements a data series called ImplicitSeries which is used by
  5. ``Plot`` class to plot implicit plots for different backends. The module,
  6. by default, implements plotting using interval arithmetic. It switches to a
  7. fall back algorithm if the expression cannot be plotted using interval arithmetic.
  8. It is also possible to specify to use the fall back algorithm for all plots.
  9. Boolean combinations of expressions cannot be plotted by the fall back
  10. algorithm.
  11. See Also
  12. ========
  13. sympy.plotting.plot
  14. References
  15. ==========
  16. .. [1] Jeffrey Allen Tupper. Reliable Two-Dimensional Graphing Methods for
  17. Mathematical Formulae with Two Free Variables.
  18. .. [2] Jeffrey Allen Tupper. Graphing Equations with Generalized Interval
  19. Arithmetic. Master's thesis. University of Toronto, 1996
  20. """
  21. from .plot import BaseSeries, Plot
  22. from .experimental_lambdify import experimental_lambdify, vectorized_lambdify
  23. from .intervalmath import interval
  24. from sympy.core.relational import (Equality, GreaterThan, LessThan,
  25. Relational, StrictLessThan, StrictGreaterThan)
  26. from sympy.core.containers import Tuple
  27. from sympy.core.relational import Eq
  28. from sympy.core.symbol import (Dummy, Symbol)
  29. from sympy.core.sympify import sympify
  30. from sympy.external import import_module
  31. from sympy.logic.boolalg import BooleanFunction
  32. from sympy.polys.polyutils import _sort_gens
  33. from sympy.utilities.decorator import doctest_depends_on
  34. from sympy.utilities.iterables import flatten
  35. import warnings
  36. class ImplicitSeries(BaseSeries):
  37. """ Representation for Implicit plot """
  38. is_implicit = True
  39. def __init__(self, expr, var_start_end_x, var_start_end_y,
  40. has_equality, use_interval_math, depth, nb_of_points,
  41. line_color):
  42. super().__init__()
  43. self.expr = sympify(expr)
  44. self.var_x = sympify(var_start_end_x[0])
  45. self.start_x = float(var_start_end_x[1])
  46. self.end_x = float(var_start_end_x[2])
  47. self.var_y = sympify(var_start_end_y[0])
  48. self.start_y = float(var_start_end_y[1])
  49. self.end_y = float(var_start_end_y[2])
  50. self.get_points = self.get_raster
  51. self.has_equality = has_equality # If the expression has equality, i.e.
  52. #Eq, Greaterthan, LessThan.
  53. self.nb_of_points = nb_of_points
  54. self.use_interval_math = use_interval_math
  55. self.depth = 4 + depth
  56. self.line_color = line_color
  57. def __str__(self):
  58. return ('Implicit equation: %s for '
  59. '%s over %s and %s over %s') % (
  60. str(self.expr),
  61. str(self.var_x),
  62. str((self.start_x, self.end_x)),
  63. str(self.var_y),
  64. str((self.start_y, self.end_y)))
  65. def get_raster(self):
  66. func = experimental_lambdify((self.var_x, self.var_y), self.expr,
  67. use_interval=True)
  68. xinterval = interval(self.start_x, self.end_x)
  69. yinterval = interval(self.start_y, self.end_y)
  70. try:
  71. func(xinterval, yinterval)
  72. except AttributeError:
  73. # XXX: AttributeError("'list' object has no attribute 'is_real'")
  74. # That needs fixing somehow - we shouldn't be catching
  75. # AttributeError here.
  76. if self.use_interval_math:
  77. warnings.warn("Adaptive meshing could not be applied to the"
  78. " expression. Using uniform meshing.", stacklevel=7)
  79. self.use_interval_math = False
  80. if self.use_interval_math:
  81. return self._get_raster_interval(func)
  82. else:
  83. return self._get_meshes_grid()
  84. def _get_raster_interval(self, func):
  85. """ Uses interval math to adaptively mesh and obtain the plot"""
  86. k = self.depth
  87. interval_list = []
  88. #Create initial 32 divisions
  89. np = import_module('numpy')
  90. xsample = np.linspace(self.start_x, self.end_x, 33)
  91. ysample = np.linspace(self.start_y, self.end_y, 33)
  92. #Add a small jitter so that there are no false positives for equality.
  93. # Ex: y==x becomes True for x interval(1, 2) and y interval(1, 2)
  94. #which will draw a rectangle.
  95. jitterx = (np.random.rand(
  96. len(xsample)) * 2 - 1) * (self.end_x - self.start_x) / 2**20
  97. jittery = (np.random.rand(
  98. len(ysample)) * 2 - 1) * (self.end_y - self.start_y) / 2**20
  99. xsample += jitterx
  100. ysample += jittery
  101. xinter = [interval(x1, x2) for x1, x2 in zip(xsample[:-1],
  102. xsample[1:])]
  103. yinter = [interval(y1, y2) for y1, y2 in zip(ysample[:-1],
  104. ysample[1:])]
  105. interval_list = [[x, y] for x in xinter for y in yinter]
  106. plot_list = []
  107. #recursive call refinepixels which subdivides the intervals which are
  108. #neither True nor False according to the expression.
  109. def refine_pixels(interval_list):
  110. """ Evaluates the intervals and subdivides the interval if the
  111. expression is partially satisfied."""
  112. temp_interval_list = []
  113. plot_list = []
  114. for intervals in interval_list:
  115. #Convert the array indices to x and y values
  116. intervalx = intervals[0]
  117. intervaly = intervals[1]
  118. func_eval = func(intervalx, intervaly)
  119. #The expression is valid in the interval. Change the contour
  120. #array values to 1.
  121. if func_eval[1] is False or func_eval[0] is False:
  122. pass
  123. elif func_eval == (True, True):
  124. plot_list.append([intervalx, intervaly])
  125. elif func_eval[1] is None or func_eval[0] is None:
  126. #Subdivide
  127. avgx = intervalx.mid
  128. avgy = intervaly.mid
  129. a = interval(intervalx.start, avgx)
  130. b = interval(avgx, intervalx.end)
  131. c = interval(intervaly.start, avgy)
  132. d = interval(avgy, intervaly.end)
  133. temp_interval_list.append([a, c])
  134. temp_interval_list.append([a, d])
  135. temp_interval_list.append([b, c])
  136. temp_interval_list.append([b, d])
  137. return temp_interval_list, plot_list
  138. while k >= 0 and len(interval_list):
  139. interval_list, plot_list_temp = refine_pixels(interval_list)
  140. plot_list.extend(plot_list_temp)
  141. k = k - 1
  142. #Check whether the expression represents an equality
  143. #If it represents an equality, then none of the intervals
  144. #would have satisfied the expression due to floating point
  145. #differences. Add all the undecided values to the plot.
  146. if self.has_equality:
  147. for intervals in interval_list:
  148. intervalx = intervals[0]
  149. intervaly = intervals[1]
  150. func_eval = func(intervalx, intervaly)
  151. if func_eval[1] and func_eval[0] is not False:
  152. plot_list.append([intervalx, intervaly])
  153. return plot_list, 'fill'
  154. def _get_meshes_grid(self):
  155. """Generates the mesh for generating a contour.
  156. In the case of equality, ``contour`` function of matplotlib can
  157. be used. In other cases, matplotlib's ``contourf`` is used.
  158. """
  159. equal = False
  160. if isinstance(self.expr, Equality):
  161. expr = self.expr.lhs - self.expr.rhs
  162. equal = True
  163. elif isinstance(self.expr, (GreaterThan, StrictGreaterThan)):
  164. expr = self.expr.lhs - self.expr.rhs
  165. elif isinstance(self.expr, (LessThan, StrictLessThan)):
  166. expr = self.expr.rhs - self.expr.lhs
  167. else:
  168. raise NotImplementedError("The expression is not supported for "
  169. "plotting in uniform meshed plot.")
  170. np = import_module('numpy')
  171. xarray = np.linspace(self.start_x, self.end_x, self.nb_of_points)
  172. yarray = np.linspace(self.start_y, self.end_y, self.nb_of_points)
  173. x_grid, y_grid = np.meshgrid(xarray, yarray)
  174. func = vectorized_lambdify((self.var_x, self.var_y), expr)
  175. z_grid = func(x_grid, y_grid)
  176. z_grid[np.ma.where(z_grid < 0)] = -1
  177. z_grid[np.ma.where(z_grid > 0)] = 1
  178. if equal:
  179. return xarray, yarray, z_grid, 'contour'
  180. else:
  181. return xarray, yarray, z_grid, 'contourf'
  182. @doctest_depends_on(modules=('matplotlib',))
  183. def plot_implicit(expr, x_var=None, y_var=None, adaptive=True, depth=0,
  184. points=300, line_color="blue", show=True, **kwargs):
  185. """A plot function to plot implicit equations / inequalities.
  186. Arguments
  187. =========
  188. - ``expr`` : The equation / inequality that is to be plotted.
  189. - ``x_var`` (optional) : symbol to plot on x-axis or tuple giving symbol
  190. and range as ``(symbol, xmin, xmax)``
  191. - ``y_var`` (optional) : symbol to plot on y-axis or tuple giving symbol
  192. and range as ``(symbol, ymin, ymax)``
  193. If neither ``x_var`` nor ``y_var`` are given then the free symbols in the
  194. expression will be assigned in the order they are sorted.
  195. The following keyword arguments can also be used:
  196. - ``adaptive`` Boolean. The default value is set to True. It has to be
  197. set to False if you want to use a mesh grid.
  198. - ``depth`` integer. The depth of recursion for adaptive mesh grid.
  199. Default value is 0. Takes value in the range (0, 4).
  200. - ``points`` integer. The number of points if adaptive mesh grid is not
  201. used. Default value is 300.
  202. - ``show`` Boolean. Default value is True. If set to False, the plot will
  203. not be shown. See ``Plot`` for further information.
  204. - ``title`` string. The title for the plot.
  205. - ``xlabel`` string. The label for the x-axis
  206. - ``ylabel`` string. The label for the y-axis
  207. Aesthetics options:
  208. - ``line_color``: float or string. Specifies the color for the plot.
  209. See ``Plot`` to see how to set color for the plots.
  210. Default value is "Blue"
  211. plot_implicit, by default, uses interval arithmetic to plot functions. If
  212. the expression cannot be plotted using interval arithmetic, it defaults to
  213. a generating a contour using a mesh grid of fixed number of points. By
  214. setting adaptive to False, you can force plot_implicit to use the mesh
  215. grid. The mesh grid method can be effective when adaptive plotting using
  216. interval arithmetic, fails to plot with small line width.
  217. Examples
  218. ========
  219. Plot expressions:
  220. .. plot::
  221. :context: reset
  222. :format: doctest
  223. :include-source: True
  224. >>> from sympy import plot_implicit, symbols, Eq, And
  225. >>> x, y = symbols('x y')
  226. Without any ranges for the symbols in the expression:
  227. .. plot::
  228. :context: close-figs
  229. :format: doctest
  230. :include-source: True
  231. >>> p1 = plot_implicit(Eq(x**2 + y**2, 5))
  232. With the range for the symbols:
  233. .. plot::
  234. :context: close-figs
  235. :format: doctest
  236. :include-source: True
  237. >>> p2 = plot_implicit(
  238. ... Eq(x**2 + y**2, 3), (x, -3, 3), (y, -3, 3))
  239. With depth of recursion as argument:
  240. .. plot::
  241. :context: close-figs
  242. :format: doctest
  243. :include-source: True
  244. >>> p3 = plot_implicit(
  245. ... Eq(x**2 + y**2, 5), (x, -4, 4), (y, -4, 4), depth = 2)
  246. Using mesh grid and not using adaptive meshing:
  247. .. plot::
  248. :context: close-figs
  249. :format: doctest
  250. :include-source: True
  251. >>> p4 = plot_implicit(
  252. ... Eq(x**2 + y**2, 5), (x, -5, 5), (y, -2, 2),
  253. ... adaptive=False)
  254. Using mesh grid without using adaptive meshing with number of points
  255. specified:
  256. .. plot::
  257. :context: close-figs
  258. :format: doctest
  259. :include-source: True
  260. >>> p5 = plot_implicit(
  261. ... Eq(x**2 + y**2, 5), (x, -5, 5), (y, -2, 2),
  262. ... adaptive=False, points=400)
  263. Plotting regions:
  264. .. plot::
  265. :context: close-figs
  266. :format: doctest
  267. :include-source: True
  268. >>> p6 = plot_implicit(y > x**2)
  269. Plotting Using boolean conjunctions:
  270. .. plot::
  271. :context: close-figs
  272. :format: doctest
  273. :include-source: True
  274. >>> p7 = plot_implicit(And(y > x, y > -x))
  275. When plotting an expression with a single variable (y - 1, for example),
  276. specify the x or the y variable explicitly:
  277. .. plot::
  278. :context: close-figs
  279. :format: doctest
  280. :include-source: True
  281. >>> p8 = plot_implicit(y - 1, y_var=y)
  282. >>> p9 = plot_implicit(x - 1, x_var=x)
  283. """
  284. has_equality = False # Represents whether the expression contains an Equality,
  285. #GreaterThan or LessThan
  286. def arg_expand(bool_expr):
  287. """
  288. Recursively expands the arguments of an Boolean Function
  289. """
  290. for arg in bool_expr.args:
  291. if isinstance(arg, BooleanFunction):
  292. arg_expand(arg)
  293. elif isinstance(arg, Relational):
  294. arg_list.append(arg)
  295. arg_list = []
  296. if isinstance(expr, BooleanFunction):
  297. arg_expand(expr)
  298. #Check whether there is an equality in the expression provided.
  299. if any(isinstance(e, (Equality, GreaterThan, LessThan))
  300. for e in arg_list):
  301. has_equality = True
  302. elif not isinstance(expr, Relational):
  303. expr = Eq(expr, 0)
  304. has_equality = True
  305. elif isinstance(expr, (Equality, GreaterThan, LessThan)):
  306. has_equality = True
  307. xyvar = [i for i in (x_var, y_var) if i is not None]
  308. free_symbols = expr.free_symbols
  309. range_symbols = Tuple(*flatten(xyvar)).free_symbols
  310. undeclared = free_symbols - range_symbols
  311. if len(free_symbols & range_symbols) > 2:
  312. raise NotImplementedError("Implicit plotting is not implemented for "
  313. "more than 2 variables")
  314. #Create default ranges if the range is not provided.
  315. default_range = Tuple(-5, 5)
  316. def _range_tuple(s):
  317. if isinstance(s, Symbol):
  318. return Tuple(s) + default_range
  319. if len(s) == 3:
  320. return Tuple(*s)
  321. raise ValueError('symbol or `(symbol, min, max)` expected but got %s' % s)
  322. if len(xyvar) == 0:
  323. xyvar = list(_sort_gens(free_symbols))
  324. var_start_end_x = _range_tuple(xyvar[0])
  325. x = var_start_end_x[0]
  326. if len(xyvar) != 2:
  327. if x in undeclared or not undeclared:
  328. xyvar.append(Dummy('f(%s)' % x.name))
  329. else:
  330. xyvar.append(undeclared.pop())
  331. var_start_end_y = _range_tuple(xyvar[1])
  332. #Check whether the depth is greater than 4 or less than 0.
  333. if depth > 4:
  334. depth = 4
  335. elif depth < 0:
  336. depth = 0
  337. series_argument = ImplicitSeries(expr, var_start_end_x, var_start_end_y,
  338. has_equality, adaptive, depth,
  339. points, line_color)
  340. #set the x and y limits
  341. kwargs['xlim'] = tuple(float(x) for x in var_start_end_x[1:])
  342. kwargs['ylim'] = tuple(float(y) for y in var_start_end_y[1:])
  343. # set the x and y labels
  344. kwargs.setdefault('xlabel', var_start_end_x[0].name)
  345. kwargs.setdefault('ylabel', var_start_end_y[0].name)
  346. p = Plot(series_argument, **kwargs)
  347. if show:
  348. p.show()
  349. return p