1# -*- coding: utf-8 -*- 2 3""" 4Plotting 5""" 6 7 8from math import sin, cos, pi, sqrt, isnan, isinf 9import numbers 10import itertools 11import palettable 12 13from mathics.version import __version__ # noqa used in loading to check consistency. 14from mathics.core.expression import ( 15 Expression, 16 Real, 17 MachineReal, 18 Symbol, 19 String, 20 Integer, 21 Integer0, 22 from_python, 23 SymbolList, 24 SymbolN, 25 SymbolRule, 26) 27 28from mathics.builtin.base import Builtin 29from mathics.builtin.graphics import Graphics 30from mathics.builtin.drawing.graphics3d import Graphics3D 31from mathics.builtin.numeric import chop 32from mathics.builtin.options import options_to_rules 33from mathics.builtin.scoping import dynamic_scoping 34 35 36try: 37 from mathics.builtin.compile import _compile, CompileArg, CompileError, real_type 38 39 has_compile = True 40except ImportError: 41 has_compile = False 42 43def gradient_palette(color_function, n, evaluation): # always returns RGB values 44 if isinstance(color_function, String): 45 color_data = Expression("ColorData", color_function).evaluate(evaluation) 46 if not color_data.has_form("ColorDataFunction", 4): 47 return 48 name, kind, interval, blend = color_data.leaves 49 if not isinstance(kind, String) or kind.get_string_value() != "Gradients": 50 return 51 if not interval.has_form("List", 2): 52 return 53 x0, x1 = (x.round_to_float() for x in interval.leaves) 54 else: 55 blend = color_function 56 x0 = 0.0 57 x1 = 1.0 58 59 xd = x1 - x0 60 offsets = [MachineReal(x0 + float(xd * i) / float(n - 1)) for i in range(n)] 61 colors = Expression("Map", blend, Expression(SymbolList, *offsets)).evaluate( 62 evaluation 63 ) 64 if len(colors.leaves) != n: 65 return 66 67 from mathics.builtin.graphics import expression_to_color, ColorError 68 69 try: 70 objects = [expression_to_color(x) for x in colors.leaves] 71 if any(x is None for x in objects): 72 return None 73 return [x.to_rgba()[:3] for x in objects] 74 except ColorError: 75 return 76 77 78class ColorDataFunction(Builtin): 79 pass 80 81 82class _GradientColorScheme(object): 83 def color_data_function(self, name): 84 colors = Expression( 85 "List", *[Expression("RGBColor", *color) for color in self.colors()] 86 ) 87 blend = Expression( 88 "Function", Expression("Blend", colors, Expression("Slot", 1)) 89 ) 90 arguments = [ 91 String(name), 92 String("Gradients"), 93 Expression(SymbolList, 0, 1), 94 blend, 95 ] 96 return Expression("ColorDataFunction", *arguments) 97 98 99class _PredefinedGradient(_GradientColorScheme): 100 def __init__(self, colors): 101 self._colors = colors 102 103 def colors(self): 104 return self._colors 105 106 107class _PalettableGradient(_GradientColorScheme): 108 def __init__(self, palette, reversed): 109 self.palette = palette 110 self.reversed = reversed 111 112 def colors(self): 113 colors = self.palette.mpl_colors 114 if self.reversed: 115 colors = list(reversed(colors)) 116 return colors 117 118 119class ColorData(Builtin): 120 """ 121 <dl> 122 <dt>'ColorData["$name$"]' 123 <dd>returns a color function with the given $name$. 124 </dl> 125 126 Define a user-defined color function: 127 >> Unprotect[ColorData]; ColorData["test"] := ColorDataFunction["test", "Gradients", {0, 1}, Blend[{Red, Green, Blue}, #1] &]; Protect[ColorData] 128 129 Compare it to the default color function, 'LakeColors': 130 >> {DensityPlot[x + y, {x, -1, 1}, {y, -1, 1}], DensityPlot[x + y, {x, -1, 1}, {y, -1, 1}, ColorFunction->"test"]} 131 = {-Graphics-, -Graphics-} 132 """ 133 134 # rules = { 135 # 'ColorData["LakeColors"]': ( 136 # """ColorDataFunction["LakeColors", "Gradients", {0, 1}, 137 # Blend[{RGBColor[0.293416, 0.0574044, 0.529412], 138 # RGBColor[0.563821, 0.527565, 0.909499], 139 # RGBColor[0.762631, 0.846998, 0.914031], 140 # RGBColor[0.941176, 0.906538, 0.834043]}, #1] &]"""), 141 # } 142 143 messages = { 144 "notent": "`1` is not a known color scheme. ColorData[] gives you lists of schemes.", 145 } 146 147 palettes = { 148 "LakeColors": _PredefinedGradient( 149 [ 150 (0.293416, 0.0574044, 0.529412), 151 (0.563821, 0.527565, 0.909499), 152 (0.762631, 0.846998, 0.914031), 153 (0.941176, 0.906538, 0.834043), 154 ] 155 ), 156 "Pastel": _PalettableGradient( 157 palettable.colorbrewer.qualitative.Pastel1_9, False 158 ), 159 "Rainbow": _PalettableGradient( 160 palettable.colorbrewer.diverging.Spectral_9, True 161 ), 162 "RedBlueTones": _PalettableGradient( 163 palettable.colorbrewer.diverging.RdBu_11, False 164 ), 165 "GreenPinkTones": _PalettableGradient( 166 palettable.colorbrewer.diverging.PiYG_9, False 167 ), 168 "GrayTones": _PalettableGradient( 169 palettable.colorbrewer.sequential.Greys_9, False 170 ), 171 "SolarColors": _PalettableGradient( 172 palettable.colorbrewer.sequential.OrRd_9, True 173 ), 174 "CherryTones": _PalettableGradient( 175 palettable.colorbrewer.sequential.Reds_9, True 176 ), 177 "FuchsiaTones": _PalettableGradient( 178 palettable.colorbrewer.sequential.RdPu_9, True 179 ), 180 "SiennaTones": _PalettableGradient( 181 palettable.colorbrewer.sequential.Oranges_9, True 182 ), 183 # specific to Mathics 184 "Paired": _PalettableGradient( 185 palettable.colorbrewer.qualitative.Paired_9, False 186 ), 187 "Accent": _PalettableGradient( 188 palettable.colorbrewer.qualitative.Accent_8, False 189 ), 190 "Aquatic": _PalettableGradient(palettable.wesanderson.Aquatic1_5, False), 191 "Zissou": _PalettableGradient(palettable.wesanderson.Zissou_5, False), 192 "Tableau": _PalettableGradient(palettable.tableau.Tableau_20, False), 193 "TrafficLight": _PalettableGradient(palettable.tableau.TrafficLight_9, False), 194 "Moonrise1": _PalettableGradient(palettable.wesanderson.Moonrise1_5, False), 195 } 196 197 def apply_directory(self, evaluation): 198 "ColorData[]" 199 return Expression(SymbolList, String("Gradients")) 200 201 def apply(self, name, evaluation): 202 "ColorData[name_String]" 203 py_name = name.get_string_value() 204 if py_name == "Gradients": 205 return Expression( 206 SymbolList, *[String(name) for name in self.palettes.keys()] 207 ) 208 palette = ColorData.palettes.get(py_name, None) 209 if palette is None: 210 evaluation.message("ColorData", "notent", name) 211 return 212 return palette.color_data_function(py_name) 213 214 @staticmethod 215 def colors(name, evaluation): 216 palette = ColorData.palettes.get(name, None) 217 if palette is None: 218 evaluation.message("ColorData", "notent", name) 219 return None 220 return palette.colors() 221 222 223def extract_pyreal(value): 224 if isinstance(value, Real): 225 return chop(value).round_to_float() 226 return None 227 228 229def zero_to_one(value): 230 if value == 0: 231 return 1 232 return value 233 234 235def compile_quiet_function(expr, arg_names, evaluation, expect_list): 236 """ 237 Given an expression return a quiet callable version. 238 Compiles the expression where possible. 239 """ 240 if has_compile and not expect_list: 241 try: 242 cfunc = _compile( 243 expr, [CompileArg(arg_name, real_type) for arg_name in arg_names] 244 ) 245 except CompileError: 246 pass 247 else: 248 249 def quiet_f(*args): 250 try: 251 result = cfunc(*args) 252 if not (isnan(result) or isinf(result)): 253 return result 254 except: 255 pass 256 return None 257 258 return quiet_f 259 260 expr = Expression(SymbolN, expr) 261 quiet_expr = Expression( 262 "Quiet", 263 expr, 264 Expression( 265 SymbolList, Expression("MessageName", Symbol("Power"), String("infy")) 266 ), 267 ) 268 269 def quiet_f(*args): 270 vars = {arg_name: Real(arg) for arg_name, arg in zip(arg_names, args)} 271 value = dynamic_scoping(quiet_expr.evaluate, vars, evaluation) 272 if expect_list: 273 if value.has_form("List", None): 274 value = [extract_pyreal(item) for item in value.leaves] 275 if any(item is None for item in value): 276 return None 277 return value 278 else: 279 return None 280 else: 281 value = extract_pyreal(value) 282 if value is None or isinf(value) or isnan(value): 283 return None 284 return value 285 286 return quiet_f 287 288 289def automatic_plot_range(values): 290 """Calculates mean and standard deviation, throwing away all points 291 which are more than 'thresh' number of standard deviations away from 292 the mean. These are then used to find good vmin and vmax values. These 293 values can then be used to find Automatic Plotrange.""" 294 295 if not values: 296 return 0, 1 297 298 thresh = 2.0 299 values = sorted(values) 300 valavg = sum(values) / len(values) 301 valdev = sqrt( 302 sum([(x - valavg) ** 2 for x in values]) / zero_to_one(len(values) - 1) 303 ) 304 305 n1, n2 = 0, len(values) - 1 306 if valdev != 0: 307 for v in values: 308 if abs(v - valavg) / valdev < thresh: 309 break 310 n1 += 1 311 for v in values[::-1]: 312 if abs(v - valavg) / valdev < thresh: 313 break 314 n2 -= 1 315 316 vrange = values[n2] - values[n1] 317 vmin = values[n1] - 0.05 * vrange # 5% extra looks nice 318 vmax = values[n2] + 0.05 * vrange 319 return vmin, vmax 320 321 322def get_plot_range(values, all_values, option): 323 if option == "System`Automatic": 324 result = automatic_plot_range(values) 325 elif option == "System`All": 326 if not all_values: 327 result = [0, 1] 328 else: 329 result = min(all_values), max(all_values) 330 else: 331 result = option 332 if result[0] == result[1]: 333 value = result[0] 334 if value > 0: 335 return 0, value * 2 336 if value < 0: 337 return value * 2, 0 338 return -1, 1 339 return result 340 341 342class _Plot(Builtin): 343 344 attributes = ("HoldAll",) 345 346 options = Graphics.options.copy() 347 options.update( 348 { 349 "Axes": "True", 350 "AspectRatio": "1 / GoldenRatio", 351 "MaxRecursion": "Automatic", 352 "Mesh": "None", 353 "PlotRange": "Automatic", 354 "PlotPoints": "None", 355 "Exclusions": "Automatic", 356 "$OptionSyntax": "Strict", 357 } 358 ) 359 360 messages = { 361 "invmaxrec": ( 362 "MaxRecursion must be a non-negative integer; the recursion value " 363 "is limited to `2`. Using MaxRecursion -> `1`." 364 ), 365 "prng": ( 366 "Value of option PlotRange -> `1` is not All, Automatic or " 367 "an appropriate list of range specifications." 368 ), 369 "ppts": "Value of option PlotPoints -> `1` is not an integer >= 2.", 370 "invexcl": ( 371 "Value of Exclusions -> `1` is not None, Automatic or an " 372 "appropriate list of constraints." 373 ), 374 } 375 376 expect_list = False 377 378 def apply(self, functions, x, start, stop, evaluation, options): 379 """%(name)s[functions_, {x_Symbol, start_, stop_}, 380 OptionsPattern[%(name)s]]""" 381 if functions.is_symbol() and functions.name is not x.get_name(): 382 rules = evaluation.definitions.get_ownvalues(functions.name) 383 for rule in rules: 384 functions = rule.apply(functions, evaluation, fully=True) 385 386 if functions.get_head_name() == "List": 387 functions_param = self.get_functions_param(functions) 388 for index, f in enumerate(functions_param): 389 if f.is_symbol() and f.name is not x.get_name(): 390 rules = evaluation.definitions.get_ownvalues(f.name) 391 for rule in rules: 392 f = rule.apply(f, evaluation, fully=True) 393 functions_param[index] = f 394 functions = functions.flatten(Symbol("List")) 395 396 expr_limits = Expression(SymbolList, x, start, stop) 397 expr = Expression( 398 self.get_name(), functions, expr_limits, *options_to_rules(options) 399 ) 400 functions = self.get_functions_param(functions) 401 x_name = x.get_name() 402 403 py_start = start.round_to_float(evaluation) 404 py_stop = stop.round_to_float(evaluation) 405 if py_start is None or py_stop is None: 406 return evaluation.message(self.get_name(), "plln", stop, expr) 407 if py_start >= py_stop: 408 return evaluation.message(self.get_name(), "plld", expr_limits) 409 start, stop = py_start, py_stop 410 411 # PlotRange Option 412 def check_range(range): 413 if range in ("System`Automatic", "System`All"): 414 return True 415 if isinstance(range, list) and len(range) == 2: 416 if isinstance(range[0], numbers.Real) and isinstance( # noqa 417 range[1], numbers.Real 418 ): 419 return True 420 return False 421 422 plotrange_option = self.get_option(options, "PlotRange", evaluation) 423 plotrange = plotrange_option.to_python(n_evaluation=evaluation) 424 x_range, y_range = self.get_plotrange(plotrange, start, stop) 425 if not check_range(x_range) or not check_range(y_range): 426 evaluation.message(self.get_name(), "prng", plotrange_option) 427 x_range, y_range = [start, stop], "Automatic" 428 429 # x_range and y_range are now either Automatic, All, or of the form [min, max] 430 assert x_range in ("System`Automatic", "System`All") or isinstance( 431 x_range, list 432 ) 433 assert y_range in ("System`Automatic", "System`All") or isinstance( 434 y_range, list 435 ) 436 437 # Mesh Option 438 mesh_option = self.get_option(options, "Mesh", evaluation) 439 mesh = mesh_option.to_python() 440 if mesh not in ["System`None", "System`Full", "System`All"]: 441 evaluation.message("Mesh", "ilevels", mesh_option) 442 mesh = "System`None" 443 444 # PlotPoints Option 445 plotpoints_option = self.get_option(options, "PlotPoints", evaluation) 446 plotpoints = plotpoints_option.to_python() 447 if plotpoints == "System`None": 448 plotpoints = 57 449 if not (isinstance(plotpoints, int) and plotpoints >= 2): 450 return evaluation.message(self.get_name(), "ppts", plotpoints) 451 452 # MaxRecursion Option 453 max_recursion_limit = 15 454 maxrecursion_option = self.get_option(options, "MaxRecursion", evaluation) 455 maxrecursion = maxrecursion_option.to_python() 456 try: 457 if maxrecursion == "System`Automatic": 458 maxrecursion = 3 459 elif maxrecursion == float("inf"): 460 maxrecursion = max_recursion_limit 461 raise ValueError 462 elif isinstance(maxrecursion, int): 463 if maxrecursion > max_recursion_limit: 464 maxrecursion = max_recursion_limit 465 raise ValueError 466 if maxrecursion < 0: 467 maxrecursion = 0 468 raise ValueError 469 else: 470 maxrecursion = 0 471 raise ValueError 472 except ValueError: 473 evaluation.message( 474 self.get_name(), "invmaxrec", maxrecursion, max_recursion_limit 475 ) 476 assert isinstance(maxrecursion, int) 477 478 # Exclusions Option 479 # TODO: Make exclusions option work properly with ParametricPlot 480 def check_exclusion(excl): 481 if isinstance(excl, list): 482 return all(check_exclusion(e) for e in excl) 483 if excl == "System`Automatic": 484 return True 485 if not isinstance(excl, numbers.Real): 486 return False 487 return True 488 489 exclusions_option = self.get_option(options, "Exclusions", evaluation) 490 exclusions = exclusions_option.to_python(n_evaluation=evaluation) 491 # TODO Turn expressions into points E.g. Sin[x] == 0 becomes 0, 2 Pi... 492 493 if exclusions in ["System`None", ["System`None"]]: 494 exclusions = "System`None" 495 elif not isinstance(exclusions, list): 496 exclusions = [exclusions] 497 498 if isinstance(exclusions, list) and all( # noqa 499 check_exclusion(excl) for excl in exclusions 500 ): 501 pass 502 503 else: 504 evaluation.message(self.get_name(), "invexcl", exclusions_option) 505 exclusions = ["System`Automatic"] 506 507 # exclusions is now either 'None' or a list of reals and 'Automatic' 508 assert exclusions == "System`None" or isinstance(exclusions, list) 509 510 # constants to generate colors 511 hue = 0.67 512 hue_pos = 0.236068 513 hue_neg = -0.763932 514 515 def get_points_minmax(points): 516 xmin = xmax = ymin = ymax = None 517 for line in points: 518 for x, y in line: 519 if xmin is None or x < xmin: 520 xmin = x 521 if xmax is None or x > xmax: 522 xmax = x 523 if ymin is None or y < ymin: 524 ymin = y 525 if ymax is None or y > ymax: 526 ymax = y 527 return xmin, xmax, ymin, ymax 528 529 def get_points_range(points): 530 xmin, xmax, ymin, ymax = get_points_minmax(points) 531 if xmin is None or xmax is None: 532 xmin, xmax = 0, 1 533 if ymin is None or ymax is None: 534 ymin, ymax = 0, 1 535 return zero_to_one(xmax - xmin), zero_to_one(ymax - ymin) 536 537 function_hues = [] 538 base_plot_points = [] # list of points in base subdivision 539 plot_points = [] # list of all plotted points 540 mesh_points = [] 541 graphics = [] # list of resulting graphics primitives 542 for index, f in enumerate(functions): 543 points = [] 544 xvalues = [] # x value for each point in points 545 tmp_mesh_points = [] # For this function only 546 continuous = False 547 d = (stop - start) / (plotpoints - 1) 548 cf = compile_quiet_function(f, [x_name], evaluation, self.expect_list) 549 for i in range(plotpoints): 550 x_value = start + i * d 551 point = self.eval_f(cf, x_value) 552 if point is not None: 553 if continuous: 554 points[-1].append(point) 555 xvalues[-1].append(x_value) 556 else: 557 points.append([point]) 558 xvalues.append([x_value]) 559 continuous = True 560 else: 561 continuous = False 562 563 base_points = [] 564 for line in points: 565 base_points.extend(line) 566 base_plot_points.extend(base_points) 567 568 xmin, xmax = automatic_plot_range([xx for xx, yy in base_points]) 569 xscale = 1.0 / zero_to_one(xmax - xmin) 570 ymin, ymax = automatic_plot_range([yy for xx, yy in base_points]) 571 yscale = 1.0 / zero_to_one(ymax - ymin) 572 573 if mesh == "System`Full": 574 for line in points: 575 tmp_mesh_points.extend(line) 576 577 def find_excl(excl): 578 # Find which line the exclusion is in 579 for l in range(len(xvalues)): # TODO: Binary Search faster? 580 if xvalues[l][0] <= excl and xvalues[l][-1] >= excl: 581 break 582 if ( 583 xvalues[l][-1] <= excl 584 and xvalues[min(l + 1, len(xvalues) - 1)][0] >= excl # nopep8 585 ): 586 return min(l + 1, len(xvalues) - 1), 0, False 587 xi = 0 588 for xi in range(len(xvalues[l]) - 1): 589 if xvalues[l][xi] <= excl and xvalues[l][xi + 1] >= excl: 590 return l, xi + 1, True 591 return l, xi + 1, False 592 593 if exclusions != "System`None": 594 for excl in exclusions: 595 if excl != "System`Automatic": 596 l, xi, split_required = find_excl(excl) 597 if split_required: 598 xvalues.insert(l + 1, xvalues[l][xi:]) 599 xvalues[l] = xvalues[l][:xi] 600 points.insert(l + 1, points[l][xi:]) 601 points[l] = points[l][:xi] 602 # assert(xvalues[l][-1] <= excl <= xvalues[l+1][0]) 603 604 # Adaptive Sampling - loop again and interpolate highly angled 605 # sections 606 607 # Cos of the maximum angle between successive line segments 608 ang_thresh = cos(pi / 180) 609 610 for line, line_xvalues in zip(points, xvalues): 611 recursion_count = 0 612 smooth = False 613 while not smooth and recursion_count < maxrecursion: 614 recursion_count += 1 615 smooth = True 616 i = 2 617 while i < len(line): 618 vec1 = ( 619 xscale * (line[i - 1][0] - line[i - 2][0]), 620 yscale * (line[i - 1][1] - line[i - 2][1]), 621 ) 622 vec2 = ( 623 xscale * (line[i][0] - line[i - 1][0]), 624 yscale * (line[i][1] - line[i - 1][1]), 625 ) 626 try: 627 angle = (vec1[0] * vec2[0] + vec1[1] * vec2[1]) / sqrt( 628 (vec1[0] ** 2 + vec1[1] ** 2) 629 * (vec2[0] ** 2 + vec2[1] ** 2) 630 ) 631 except ZeroDivisionError: 632 angle = 0.0 633 if abs(angle) < ang_thresh: 634 smooth = False 635 incr = 0 636 637 x_value = 0.5 * (line_xvalues[i - 1] + line_xvalues[i]) 638 639 point = self.eval_f(cf, x_value) 640 if point is not None: 641 line.insert(i, point) 642 line_xvalues.insert(i, x_value) 643 incr += 1 644 645 x_value = 0.5 * (line_xvalues[i - 2] + line_xvalues[i - 1]) 646 point = self.eval_f(cf, x_value) 647 if point is not None: 648 line.insert(i - 1, point) 649 line_xvalues.insert(i - 1, x_value) 650 incr += 1 651 652 i += incr 653 i += 1 654 655 if exclusions == "System`None": # Join all the Lines 656 points = [[(xx, yy) for line in points for xx, yy in line]] 657 658 graphics.append(Expression("Hue", hue, 0.6, 0.6)) 659 graphics.append(Expression("Line", from_python(points))) 660 661 for line in points: 662 plot_points.extend(line) 663 664 if mesh == "System`All": 665 for line in points: 666 tmp_mesh_points.extend(line) 667 668 if mesh != "System`None": 669 mesh_points.append(tmp_mesh_points) 670 671 function_hues.append(hue) 672 673 if index % 4 == 0: 674 hue += hue_pos 675 else: 676 hue += hue_neg 677 if hue > 1: 678 hue -= 1 679 if hue < 0: 680 hue += 1 681 682 x_range = get_plot_range( 683 [xx for xx, yy in base_plot_points], [xx for xx, yy in plot_points], x_range 684 ) 685 y_range = get_plot_range( 686 [yy for xx, yy in base_plot_points], [yy for xx, yy in plot_points], y_range 687 ) 688 689 options["System`PlotRange"] = from_python([x_range, y_range]) 690 691 if mesh != "None": 692 for hue, points in zip(function_hues, mesh_points): 693 graphics.append(Expression("Hue", hue, 0.6, 0.6)) 694 meshpoints = [Expression(SymbolList, xx, yy) for xx, yy in points] 695 graphics.append( 696 Expression("Point", Expression(SymbolList, *meshpoints)) 697 ) 698 699 return Expression( 700 "Graphics", Expression(SymbolList, *graphics), *options_to_rules(options) 701 ) 702 703 704class _Chart(Builtin): 705 attributes = ("HoldAll",) 706 707 options = Graphics.options.copy() 708 options.update( 709 { 710 "Mesh": "None", 711 "PlotRange": "Automatic", 712 "ChartLabels": "None", 713 "ChartLegends": "None", 714 "ChartStyle": "Automatic", 715 } 716 ) 717 718 never_monochrome = False 719 720 def _draw(self, data, color, evaluation, options): 721 raise NotImplementedError() 722 723 def apply(self, points, evaluation, options): 724 "%(name)s[points_, OptionsPattern[%(name)s]]" 725 726 points = points.evaluate(evaluation) 727 728 if points.get_head_name() != "System`List" or not points.leaves: 729 return 730 731 if points.leaves[0].get_head_name() == "System`List": 732 if not all( 733 group.get_head_name() == "System`List" for group in points.leaves 734 ): 735 return 736 multiple_colors = True 737 groups = points.leaves 738 else: 739 multiple_colors = False 740 groups = [points] 741 742 chart_legends = self.get_option(options, "ChartLegends", evaluation) 743 has_chart_legends = chart_legends.get_head_name() == "System`List" 744 if has_chart_legends: 745 multiple_colors = True 746 747 def to_number(x): 748 if isinstance(x, Integer): 749 return float(x.get_int_value()) 750 return x.round_to_float(evaluation=evaluation) 751 752 data = [[to_number(x) for x in group.leaves] for group in groups] 753 754 chart_style = self.get_option(options, "ChartStyle", evaluation) 755 if ( 756 isinstance(chart_style, Symbol) 757 and chart_style.get_name() == "System`Automatic" 758 ): 759 chart_style = String("Automatic") 760 761 if chart_style.get_head_name() == "System`List": 762 colors = chart_style.leaves 763 spread_colors = False 764 elif isinstance(chart_style, String): 765 if chart_style.get_string_value() == "Automatic": 766 mpl_colors = palettable.wesanderson.Moonrise1_5.mpl_colors 767 else: 768 mpl_colors = ColorData.colors(chart_style.get_string_value()) 769 if mpl_colors is None: 770 return 771 multiple_colors = True 772 773 if not multiple_colors and not self.never_monochrome: 774 colors = [Expression("RGBColor", *mpl_colors[0])] 775 else: 776 colors = [Expression("RGBColor", *c) for c in mpl_colors] 777 spread_colors = True 778 else: 779 return 780 781 def legends(names): 782 if not data: 783 return 784 785 n = len(data[0]) 786 for d in data[1:]: 787 if len(d) != n: 788 return # data groups should have same size 789 790 def box(color): 791 return Expression( 792 "Graphics", 793 Expression( 794 "List", Expression("FaceForm", color), Expression("Rectangle") 795 ), 796 Expression( 797 SymbolRule, Symbol("ImageSize"), Expression(SymbolList, 50, 50) 798 ), 799 ) 800 801 rows_per_col = 5 802 803 n_cols = 1 + len(names) // rows_per_col 804 if len(names) % rows_per_col == 0: 805 n_cols -= 1 806 807 if n_cols == 1: 808 n_rows = len(names) 809 else: 810 n_rows = rows_per_col 811 812 for i in range(n_rows): 813 items = [] 814 for j in range(n_cols): 815 k = 1 + i + j * rows_per_col 816 if k - 1 < len(names): 817 items.extend([box(color(k, n)), names[k - 1]]) 818 else: 819 items.extend([String(""), String("")]) 820 yield Expression(SymbolList, *items) 821 822 def color(k, n): 823 if spread_colors and n < len(colors): 824 index = int(k * (len(colors) - 1)) // n 825 return colors[index] 826 else: 827 return colors[(k - 1) % len(colors)] 828 829 chart = self._draw(data, color, evaluation, options) 830 831 if has_chart_legends: 832 grid = Expression( 833 "Grid", Expression(SymbolList, *list(legends(chart_legends.leaves))) 834 ) 835 chart = Expression("Row", Expression(SymbolList, chart, grid)) 836 837 return chart 838 839 840class PieChart(_Chart): 841 """ 842 <dl> 843 <dt>'PieChart[{$p1$, $p2$ ...}]' 844 <dd>draws a pie chart. 845 </dl> 846 847 >> PieChart[{1, 4, 2}] 848 = -Graphics- 849 850 >> PieChart[{8, 16, 2}, SectorOrigin -> {Automatic, 1.5}] 851 = -Graphics- 852 853 >> PieChart[{{10, 20, 30}, {15, 22, 30}}] 854 = -Graphics- 855 856 >> PieChart[{{10, 20, 30}, {15, 22, 30}}, SectorSpacing -> None] 857 = -Graphics- 858 859 >> PieChart[{{10, 20, 30}, {15, 22, 30}}, ChartLabels -> {a, b, c}] 860 = -Graphics- 861 862 Negative values are clipped to 0. 863 >> PieChart[{1, -1, 3}] 864 = -Graphics- 865 """ 866 867 options = _Chart.options.copy() 868 options.update( 869 { 870 "Axes": "{False, False}", 871 "AspectRatio": "1", 872 "SectorOrigin": "{Automatic, 0}", 873 "SectorSpacing": "Automatic", 874 } 875 ) 876 877 never_monochrome = True 878 879 def _draw(self, data, color, evaluation, options): 880 data = [[max(0.0, x) for x in group] for group in data] 881 882 sector_origin = self.get_option(options, "SectorOrigin", evaluation) 883 if not sector_origin.has_form("List", 2): 884 return 885 sector_origin = Expression(SymbolN, sector_origin).evaluate(evaluation) 886 887 orientation = sector_origin.leaves[0] 888 if ( 889 isinstance(orientation, Symbol) 890 and orientation.get_name() == "System`Automatic" 891 ): 892 sector_phi = pi 893 sector_sign = -1.0 894 elif orientation.has_form("List", 2) and isinstance( 895 orientation.leaves[1], String 896 ): 897 sector_phi = orientation.leaves[0].round_to_float() 898 clock_name = orientation.leaves[1].get_string_value() 899 if clock_name == "Clockwise": 900 sector_sign = -1.0 901 elif clock_name == "Counterclockwise": 902 sector_sign = 1.0 903 else: 904 return 905 else: 906 return 907 908 sector_spacing = self.get_option(options, "SectorSpacing", evaluation) 909 if isinstance(sector_spacing, Symbol): 910 if sector_spacing.get_name() == "System`Automatic": 911 sector_spacing = Expression(SymbolList, Integer0, Real(0.2)) 912 elif sector_spacing.get_name() == "System`None": 913 sector_spacing = Expression(SymbolList, Integer0, Integer0) 914 else: 915 return 916 if not sector_spacing.has_form("List", 2): 917 return 918 segment_spacing = 0.0 # not yet implemented; needs real arc graphics 919 radius_spacing = max(0.0, min(1.0, sector_spacing.leaves[1].round_to_float())) 920 921 def vector2(x, y): 922 return Expression(SymbolList, Real(x), Real(y)) 923 924 def radii(): 925 outer = 2.0 926 inner = sector_origin.leaves[1].round_to_float() 927 n = len(data) 928 929 d = (outer - inner) / n 930 931 r0 = outer 932 for i in range(n): 933 r1 = r0 - d 934 if i > 0: 935 r0 -= radius_spacing * d 936 yield (r0, r1) 937 r0 = r1 938 939 def phis(values): 940 s = sum(values) 941 942 t = 0.0 943 pi2 = pi * 2.0 944 phi0 = pi 945 spacing = sector_sign * segment_spacing / 2.0 946 947 for k, value in enumerate(values): 948 t += value 949 phi1 = sector_phi + sector_sign * (t / s) * pi2 950 951 yield (phi0 + spacing, phi1 - spacing) 952 phi0 = phi1 953 954 def segments(): 955 yield Expression("EdgeForm", Symbol("Black")) 956 957 origin = vector2(0.0, 0.0) 958 959 for values, (r0, r1) in zip(data, radii()): 960 radius = vector2(r0, r0) 961 962 n = len(values) 963 964 for k, (phi0, phi1) in enumerate(phis(values)): 965 yield Expression( 966 "Style", 967 Expression("Disk", origin, radius, vector2(phi0, phi1)), 968 color(k + 1, n), 969 ) 970 971 if r1 > 0.0: 972 yield Expression( 973 "Style", 974 Expression("Disk", origin, vector2(r1, r1)), 975 Symbol("White"), 976 ) 977 978 def labels(names): 979 yield Expression("FaceForm", Symbol("Black")) 980 981 for values, (r0, r1) in zip(data, radii()): 982 for name, (phi0, phi1) in zip(names, phis(values)): 983 r = (r0 + r1) / 2.0 984 phi = (phi0 + phi1) / 2.0 985 yield Expression("Text", name, vector2(r * cos(phi), r * sin(phi))) 986 987 graphics = list(segments()) 988 989 chart_labels = self.get_option(options, "ChartLabels", evaluation) 990 if chart_labels.get_head_name() == "System`List": 991 graphics.extend(list(labels(chart_labels.leaves))) 992 993 options["System`PlotRange"] = Expression( 994 "List", vector2(-2.0, 2.0), vector2(-2.0, 2.0) 995 ) 996 997 return Expression( 998 "Graphics", Expression(SymbolList, *graphics), *options_to_rules(options) 999 ) 1000 1001 1002class BarChart(_Chart): 1003 """ 1004 <dl> 1005 <dt>'BarChart[{$b1$, $b2$ ...}]' 1006 <dd>makes a bar chart with lengths $b1$, $b2$, .... 1007 </dl> 1008 1009 >> BarChart[{1, 4, 2}] 1010 = -Graphics- 1011 1012 >> BarChart[{1, 4, 2}, ChartStyle -> {Red, Green, Blue}] 1013 = -Graphics- 1014 1015 >> BarChart[{{1, 2, 3}, {2, 3, 4}}] 1016 = -Graphics- 1017 1018 >> BarChart[{{1, 2, 3}, {2, 3, 4}}, ChartLabels -> {"a", "b", "c"}] 1019 = -Graphics- 1020 1021 >> BarChart[{{1, 5}, {3, 4}}, ChartStyle -> {{EdgeForm[Thin], White}, {EdgeForm[Thick], White}}] 1022 = -Graphics- 1023 """ 1024 1025 options = _Chart.options.copy() 1026 options.update( 1027 {"Axes": "{False, True}", "AspectRatio": "1 / GoldenRatio",} 1028 ) 1029 1030 def _draw(self, data, color, evaluation, options): 1031 def vector2(x, y): 1032 return Expression(SymbolList, Real(x), Real(y)) 1033 1034 def boxes(): 1035 w = 0.9 1036 s = 0.06 1037 w_half = 0.5 * w 1038 x = 0.1 + s + w_half 1039 1040 for y_values in data: 1041 y_length = len(y_values) 1042 for i, y in enumerate(y_values): 1043 x0 = x - w_half 1044 x1 = x0 + w 1045 yield (i + 1, y_length), x0, x1, y 1046 x = x1 + s + w_half 1047 1048 x += 0.2 1049 1050 def rectangles(): 1051 yield Expression("EdgeForm", Symbol("Black")) 1052 1053 last_x1 = 0 1054 1055 for (k, n), x0, x1, y in boxes(): 1056 yield Expression( 1057 "Style", 1058 Expression( 1059 "Rectangle", 1060 Expression(SymbolList, x0, 0), 1061 Expression(SymbolList, x1, y), 1062 ), 1063 color(k, n), 1064 ) 1065 1066 last_x1 = x1 1067 1068 yield Expression( 1069 "Line", Expression(SymbolList, vector2(0, 0), vector2(last_x1, 0)) 1070 ) 1071 1072 def axes(): 1073 yield Expression("FaceForm", Symbol("Black")) 1074 1075 def points(x): 1076 return Expression(SymbolList, vector2(x, 0), vector2(x, -0.2)) 1077 1078 for (k, n), x0, x1, y in boxes(): 1079 if k == 1: 1080 yield Expression("Line", points(x0)) 1081 if k == n: 1082 yield Expression("Line", points(x1)) 1083 1084 def labels(names): 1085 yield Expression("FaceForm", Symbol("Black")) 1086 1087 for (k, n), x0, x1, y in boxes(): 1088 if k <= len(names): 1089 name = names[k - 1] 1090 yield Expression("Text", name, vector2((x0 + x1) / 2, -0.2)) 1091 1092 x_coords = list(itertools.chain(*[[x0, x1] for (k, n), x0, x1, y in boxes()])) 1093 y_coords = [0] + [y for (k, n), x0, x1, y in boxes()] 1094 1095 graphics = list(rectangles()) + list(axes()) 1096 1097 x_range = "System`All" 1098 y_range = "System`All" 1099 1100 x_range = list(get_plot_range(x_coords, x_coords, x_range)) 1101 y_range = list(get_plot_range(y_coords, y_coords, y_range)) 1102 1103 chart_labels = self.get_option(options, "ChartLabels", evaluation) 1104 if chart_labels.get_head_name() == "System`List": 1105 graphics.extend(list(labels(chart_labels.leaves))) 1106 y_range[0] = -0.4 # room for labels at the bottom 1107 1108 # always specify -.1 as the minimum x plot range, as this will make the y axis apppear 1109 # at origin (0,0); otherwise it will be shifted right; see GraphicsBox.axis_ticks(). 1110 x_range[0] = -0.1 1111 1112 options["System`PlotRange"] = Expression( 1113 "List", vector2(*x_range), vector2(*y_range) 1114 ) 1115 1116 return Expression( 1117 "Graphics", Expression(SymbolList, *graphics), *options_to_rules(options) 1118 ) 1119 1120 1121class Histogram(Builtin): 1122 """ 1123 <dl> 1124 <dt>'Histogram[{$x1$, $x2$ ...}]' 1125 <dd>plots a histogram using the values $x1$, $x2$, .... 1126 </dl> 1127 1128 >> Histogram[{3, 8, 10, 100, 1000, 500, 300, 200, 10, 20, 200, 100, 200, 300, 500}] 1129 = -Graphics- 1130 1131 >> Histogram[{{1, 2, 10, 5, 50, 20}, {90, 100, 101, 120, 80}}] 1132 = -Graphics- 1133 """ 1134 1135 attributes = ("HoldAll",) 1136 1137 options = Graphics.options.copy() 1138 options.update( 1139 { 1140 "Axes": "{True, True}", 1141 "AspectRatio": "1 / GoldenRatio", 1142 "Mesh": "None", 1143 "PlotRange": "Automatic", 1144 } 1145 ) 1146 1147 def apply(self, points, spec, evaluation, options): 1148 "%(name)s[points_, spec___, OptionsPattern[%(name)s]]" 1149 1150 points = points.evaluate(evaluation) 1151 spec = spec.evaluate(evaluation).get_sequence() 1152 1153 if spec and len(spec) not in (1, 2): 1154 return 1155 1156 if points.get_head_name() != "System`List" or not points.leaves: 1157 return 1158 1159 if points.leaves[0].get_head_name() == "System`List": 1160 if not all(q.get_head_name() == "System`List" for q in points.leaves): 1161 return 1162 input = points.leaves 1163 else: 1164 input = [points] 1165 1166 def to_numbers(l): 1167 for x in l: 1168 y = x.to_mpmath() 1169 if y is not None: 1170 yield y 1171 1172 matrix = [list(to_numbers(data.leaves)) for data in input] 1173 minima = [min(data) for data in matrix] 1174 maxima = [max(data) for data in matrix] 1175 max_bins = max(len(data) for data in matrix) 1176 1177 minimum = min(minima) 1178 maximum = max(maxima) 1179 1180 if minimum > 0: 1181 minimum = 0 1182 1183 span = maximum - minimum 1184 1185 from math import ceil 1186 from mpmath import floor as mpfloor, ceil as mpceil 1187 1188 class Distribution: 1189 def __init__(self, data, n_bins): 1190 bin_width = span / n_bins 1191 bins = [0] * n_bins 1192 for x in data: 1193 b = int(mpfloor((x - minimum) / bin_width)) 1194 if b < 0: 1195 b = 0 1196 elif b >= n_bins: 1197 b = n_bins - 1 1198 bins[b] += 1 1199 self.bins = bins 1200 self.bin_width = bin_width 1201 1202 def n_bins(self): 1203 return len(self.bins) 1204 1205 def cost(self): 1206 # see http://toyoizumilab.brain.riken.jp/hideaki/res/histogram.html 1207 bins = self.bins 1208 n_bins = len(bins) 1209 k = sum(bins) / n_bins 1210 v = sum(x * x for x in ((b - k) for b in bins)) / n_bins 1211 bin_width = self.bin_width 1212 return (2 * k - v) / (bin_width * bin_width) 1213 1214 def graphics(self, color): 1215 bins = self.bins 1216 n_bins = len(bins) 1217 bin_width = self.bin_width 1218 1219 def boxes(): 1220 x = minimum 1221 1222 for i, count in enumerate(bins): 1223 x1 = x + bin_width 1224 yield x, x1, count 1225 x = minimum + ((i + 1) * span) / n_bins 1226 1227 def rectangles(): 1228 yield Expression("EdgeForm", Expression("RGBColor", 0, 0, 0)) 1229 1230 last_x1 = 0 1231 style = Expression("RGBColor", *color) 1232 1233 for x0, x1, y in boxes(): 1234 yield Expression( 1235 "Style", 1236 Expression( 1237 "Rectangle", 1238 Expression(SymbolList, x0, 0), 1239 Expression(SymbolList, x1, y), 1240 ), 1241 style, 1242 ) 1243 1244 last_x1 = x1 1245 1246 yield Expression( 1247 "Line", 1248 Expression( 1249 "List", 1250 Expression(SymbolList, 0, 0), 1251 Expression(SymbolList, last_x1, 0), 1252 ), 1253 ) 1254 1255 return list(rectangles()) 1256 1257 def compute_cost(n_bins): 1258 distributions = [Distribution(data, n_bins) for data in matrix] 1259 return sum(d.cost() for d in distributions), distributions 1260 1261 def best_distributions(n_bins, dir, cost0, distributions0): 1262 if dir > 0: 1263 step_size = (max_bins - n_bins) // 2 1264 else: 1265 step_size = (n_bins - 1) // 2 1266 if step_size < 1: 1267 step_size = 1 1268 1269 while True: 1270 new_n_bins = n_bins + dir * step_size 1271 if new_n_bins < 1 or new_n_bins > max_bins: 1272 good = False 1273 else: 1274 cost, distributions = compute_cost(new_n_bins) 1275 good = cost < cost0 1276 1277 if not good: 1278 if step_size == 1: 1279 break 1280 step_size = max(step_size // 2, 1) 1281 else: 1282 n_bins = new_n_bins 1283 cost0 = cost 1284 distributions0 = distributions 1285 1286 return cost0, distributions0 1287 1288 def graphics(distributions): 1289 palette = palettable.wesanderson.FantasticFox1_5 1290 colors = list(reversed(palette.mpl_colors)) 1291 1292 from itertools import chain 1293 1294 n_bins = distributions[0].n_bins() 1295 x_coords = [minimum + (i * span) / n_bins for i in range(n_bins + 1)] 1296 y_coords = [0] + list( 1297 chain(*[distribution.bins for distribution in distributions]) 1298 ) 1299 1300 graphics = [] 1301 for i, distribution in enumerate(distributions): 1302 color = colors[i % len(colors)] 1303 graphics.extend(list(chain(*[distribution.graphics(color)]))) 1304 1305 x_range = "System`All" 1306 y_range = "System`All" 1307 1308 x_range = list(get_plot_range(x_coords, x_coords, x_range)) 1309 y_range = list(get_plot_range(y_coords, y_coords, y_range)) 1310 1311 # always specify -.1 as the minimum x plot range, as this will make the y axis apppear 1312 # at origin (0,0); otherwise it will be shifted right; see GraphicsBox.axis_ticks(). 1313 x_range[0] = -0.1 1314 1315 options["System`PlotRange"] = from_python([x_range, y_range]) 1316 1317 return Expression( 1318 "Graphics", 1319 Expression(SymbolList, *graphics), 1320 *options_to_rules(options) 1321 ) 1322 1323 def manual_bins(bspec, hspec): 1324 if isinstance(bspec, Integer): 1325 distributions = [ 1326 Distribution(data, bspec.get_int_value()) for data in matrix 1327 ] 1328 return graphics(distributions) 1329 elif bspec.get_head_name() == "System`List" and len(bspec.leaves) == 1: 1330 bin_width = bspec[0].to_mpmath() 1331 distributions = [ 1332 Distribution(data, int(mpceil(span / bin_width))) for data in matrix 1333 ] 1334 return graphics(distributions) 1335 1336 def auto_bins(): 1337 # start with Rice's rule, see https://en.wikipedia.org/wiki/Histogram 1338 n_bins = int(ceil(2 * (max_bins ** (1.0 / 3.0)))) 1339 1340 # now optimize the bin size by going into both directions and looking 1341 # for local minima. 1342 1343 cost0, distributions0 = compute_cost(n_bins) 1344 cost_r, distributions_r = best_distributions( 1345 n_bins, 1, cost0, distributions0 1346 ) 1347 cost_l, distributions_l = best_distributions( 1348 n_bins, -1, cost0, distributions0 1349 ) 1350 1351 if cost_r < cost_l: 1352 distributions = distributions_r 1353 else: 1354 distributions = distributions_l 1355 1356 return graphics(distributions) 1357 1358 if not spec: 1359 return auto_bins() 1360 else: 1361 if len(spec) < 2: 1362 spec.append(None) 1363 return manual_bins(*spec) 1364 return Expression( 1365 "Graphics", 1366 Expression(SymbolList, *graphics), 1367 *options_to_rules(options, Graphics.options) 1368 ) 1369 1370 1371class _ListPlot(Builtin): 1372 messages = { 1373 "prng": ( 1374 "Value of option PlotRange -> `1` is not All, Automatic or " 1375 "an appropriate list of range specifications." 1376 ), 1377 "joind": "Value of option Joined -> `1` is not True or False.", 1378 } 1379 1380 def apply(self, points, evaluation, options): 1381 "%(name)s[points_, OptionsPattern[%(name)s]]" 1382 1383 plot_name = self.get_name() 1384 all_points = points.to_python(n_evaluation=evaluation) 1385 expr = Expression(self.get_name(), points, *options_to_rules(options)) 1386 1387 # PlotRange Option 1388 def check_range(range): 1389 if range in ("System`Automatic", "System`All"): 1390 return True 1391 if isinstance(range, list) and len(range) == 2: 1392 if isinstance(range[0], numbers.Real) and isinstance( # noqa 1393 range[1], numbers.Real 1394 ): 1395 return True 1396 return False 1397 1398 plotrange_option = self.get_option(options, "PlotRange", evaluation) 1399 plotrange = plotrange_option.to_python(n_evaluation=evaluation) 1400 if plotrange == "System`All": 1401 plotrange = ["System`All", "System`All"] 1402 elif plotrange == "System`Automatic": 1403 plotrange = ["System`Automatic", "System`Automatic"] 1404 elif isinstance(plotrange, numbers.Real): 1405 plotrange = [[-plotrange, plotrange], [-plotrange, plotrange]] 1406 elif isinstance(plotrange, list) and len(plotrange) == 2: 1407 if all(isinstance(pr, numbers.Real) for pr in plotrange): 1408 plotrange = ["System`All", plotrange] 1409 elif all(check_range(pr) for pr in plotrange): 1410 pass 1411 else: 1412 evaluation.message(self.get_name(), "prng", plotrange_option) 1413 plotrange = ["System`Automatic", "System`Automatic"] 1414 1415 x_range, y_range = plotrange[0], plotrange[1] 1416 assert x_range in ("System`Automatic", "System`All") or isinstance( 1417 x_range, list 1418 ) 1419 assert y_range in ("System`Automatic", "System`All") or isinstance( 1420 y_range, list 1421 ) 1422 1423 # Filling option 1424 # TODO: Fill between corresponding points in two datasets: 1425 filling_option = self.get_option(options, "Filling", evaluation) 1426 filling = filling_option.to_python(n_evaluation=evaluation) 1427 if filling in [ 1428 "System`Top", 1429 "System`Bottom", 1430 "System`Axis", 1431 ] or isinstance( # noqa 1432 filling, numbers.Real 1433 ): 1434 pass 1435 else: 1436 # Mathematica does not even check that filling is sane 1437 filling = None 1438 1439 # Joined Option 1440 joined_option = self.get_option(options, "Joined", evaluation) 1441 joined = joined_option.to_python() 1442 if joined not in [True, False]: 1443 evaluation.message(plot_name, "joind", joined_option, expr) 1444 joined = False 1445 1446 if isinstance(all_points, list) and len(all_points) != 0: 1447 if all(not isinstance(point, list) for point in all_points): 1448 # Only y values given 1449 all_points = [ 1450 [[float(i + 1), all_points[i]] for i in range(len(all_points))] 1451 ] 1452 elif all(isinstance(line, list) and len(line) == 2 for line in all_points): 1453 # Single list of (x,y) pairs 1454 all_points = [all_points] 1455 elif all(isinstance(line, list) for line in all_points): 1456 # List of lines 1457 if all( 1458 isinstance(point, list) and len(point) == 2 1459 for line in all_points 1460 for point in line 1461 ): 1462 pass 1463 elif all( 1464 not isinstance(point, list) for line in all_points for point in line 1465 ): 1466 all_points = [ 1467 [[float(i + 1), l] for i, l in enumerate(line)] 1468 for line in all_points 1469 ] 1470 else: 1471 return 1472 else: 1473 return 1474 else: 1475 return 1476 1477 # Split into segments at missing data 1478 all_points = [[line] for line in all_points] 1479 for l, line in enumerate(all_points): 1480 i = 0 1481 while i < len(all_points[l]): 1482 seg = line[i] 1483 for j, point in enumerate(seg): 1484 if not ( 1485 isinstance(point[0], (int, float)) 1486 and isinstance(point[1], (int, float)) 1487 ): 1488 all_points[l].insert(i, seg[:j]) 1489 all_points[l][i + 1] = seg[j + 1 :] 1490 i -= 1 1491 break 1492 1493 i += 1 1494 1495 y_range = get_plot_range( 1496 [y for line in all_points for seg in line for x, y in seg], 1497 [y for line in all_points for seg in line for x, y in seg], 1498 y_range, 1499 ) 1500 x_range = get_plot_range( 1501 [x for line in all_points for seg in line for x, y in seg], 1502 [x for line in all_points for seg in line for x, y in seg], 1503 x_range, 1504 ) 1505 1506 if filling == "System`Axis": 1507 # TODO: Handle arbitary axis intercepts 1508 filling = 0.0 1509 elif filling == "System`Bottom": 1510 filling = y_range[0] 1511 elif filling == "System`Top": 1512 filling = y_range[1] 1513 1514 hue = 0.67 1515 hue_pos = 0.236068 1516 hue_neg = -0.763932 1517 1518 graphics = [] 1519 for indx, line in enumerate(all_points): 1520 graphics.append(Expression("Hue", hue, 0.6, 0.6)) 1521 for segment in line: 1522 if joined: 1523 graphics.append(Expression("Line", from_python(segment))) 1524 if filling is not None: 1525 graphics.append(Expression("Hue", hue, 0.6, 0.6, 0.2)) 1526 fill_area = list(segment) 1527 fill_area.append([segment[-1][0], filling]) 1528 fill_area.append([segment[0][0], filling]) 1529 graphics.append(Expression("Polygon", from_python(fill_area))) 1530 else: 1531 graphics.append(Expression("Point", from_python(segment))) 1532 if filling is not None: 1533 for point in segment: 1534 graphics.append( 1535 Expression( 1536 "Line", 1537 from_python( 1538 [[point[0], filling], [point[0], point[1]]] 1539 ), 1540 ) 1541 ) 1542 1543 if indx % 4 == 0: 1544 hue += hue_pos 1545 else: 1546 hue += hue_neg 1547 if hue > 1: 1548 hue -= 1 1549 if hue < 0: 1550 hue += 1 1551 1552 options["System`PlotRange"] = from_python([x_range, y_range]) 1553 1554 return Expression( 1555 "Graphics", 1556 Expression(SymbolList, *graphics), 1557 *options_to_rules(options, Graphics.options) 1558 ) 1559 1560 1561class _Plot3D(Builtin): 1562 messages = { 1563 "invmaxrec": ( 1564 "MaxRecursion must be a non-negative integer; the recursion value " 1565 "is limited to `2`. Using MaxRecursion -> `1`." 1566 ), 1567 "prng": ( 1568 "Value of option PlotRange -> `1` is not All, Automatic or " 1569 "an appropriate list of range specifications." 1570 ), 1571 "invmesh": "Mesh must be one of {None, Full, All}. Using Mesh->None.", 1572 "invpltpts": ( 1573 "Value of PlotPoints -> `1` is not a positive integer " 1574 "or appropriate list of positive integers." 1575 ), 1576 } 1577 1578 def apply(self, functions, x, xstart, xstop, y, ystart, ystop, evaluation, options): 1579 """%(name)s[functions_, {x_Symbol, xstart_, xstop_}, 1580 {y_Symbol, ystart_, ystop_}, OptionsPattern[%(name)s]]""" 1581 xexpr_limits = Expression(SymbolList, x, xstart, xstop) 1582 yexpr_limits = Expression(SymbolList, y, ystart, ystop) 1583 expr = Expression( 1584 self.get_name(), 1585 functions, 1586 xexpr_limits, 1587 yexpr_limits, 1588 *options_to_rules(options) 1589 ) 1590 1591 functions = self.get_functions_param(functions) 1592 plot_name = self.get_name() 1593 1594 def convert_limit(value, limits): 1595 result = value.round_to_float(evaluation) 1596 if result is None: 1597 evaluation.message(plot_name, "plln", value, limits) 1598 return result 1599 1600 xstart = convert_limit(xstart, xexpr_limits) 1601 xstop = convert_limit(xstop, xexpr_limits) 1602 ystart = convert_limit(ystart, yexpr_limits) 1603 ystop = convert_limit(ystop, yexpr_limits) 1604 if None in (xstart, xstop, ystart, ystop): 1605 return 1606 1607 if ystart >= ystop: 1608 evaluation.message(plot_name, "plln", ystop, expr) 1609 return 1610 1611 if xstart >= xstop: 1612 evaluation.message(plot_name, "plln", xstop, expr) 1613 return 1614 1615 # Mesh Option 1616 mesh_option = self.get_option(options, "Mesh", evaluation) 1617 mesh = mesh_option.to_python() 1618 if mesh not in ["System`None", "System`Full", "System`All"]: 1619 evaluation.message("Mesh", "ilevels", mesh_option) 1620 mesh = "System`Full" 1621 1622 # PlotPoints Option 1623 plotpoints_option = self.get_option(options, "PlotPoints", evaluation) 1624 plotpoints = plotpoints_option.to_python() 1625 1626 def check_plotpoints(steps): 1627 if isinstance(steps, int) and steps > 0: 1628 return True 1629 return False 1630 1631 if plotpoints == "System`None": 1632 plotpoints = [7, 7] 1633 elif check_plotpoints(plotpoints): 1634 plotpoints = [plotpoints, plotpoints] 1635 1636 if not ( 1637 isinstance(plotpoints, list) 1638 and len(plotpoints) == 2 1639 and check_plotpoints(plotpoints[0]) 1640 and check_plotpoints(plotpoints[1]) 1641 ): 1642 evaluation.message(self.get_name(), "invpltpts", plotpoints) 1643 plotpoints = [7, 7] 1644 1645 # MaxRecursion Option 1646 maxrec_option = self.get_option(options, "MaxRecursion", evaluation) 1647 max_depth = maxrec_option.to_python() 1648 if isinstance(max_depth, int): 1649 if max_depth < 0: 1650 max_depth = 0 1651 evaluation.message(self.get_name(), "invmaxrec", max_depth, 15) 1652 elif max_depth > 15: 1653 max_depth = 15 1654 evaluation.message(self.get_name(), "invmaxrec", max_depth, 15) 1655 else: 1656 pass # valid 1657 elif max_depth == float("inf"): 1658 max_depth = 15 1659 evaluation.message(self.get_name(), "invmaxrec", max_depth, 15) 1660 else: 1661 max_depth = 0 1662 evaluation.message(self.get_name(), "invmaxrec", max_depth, 15) 1663 1664 # Plot the functions 1665 graphics = [] 1666 for indx, f in enumerate(functions): 1667 stored = {} 1668 1669 cf = compile_quiet_function( 1670 f, [x.get_name(), y.get_name()], evaluation, False 1671 ) 1672 1673 def eval_f(x_value, y_value): 1674 try: 1675 return stored[(x_value, y_value)] 1676 except KeyError: 1677 value = cf(x_value, y_value) 1678 if value is not None: 1679 value = float(value) 1680 stored[(x_value, y_value)] = value 1681 return value 1682 1683 triangles = [] 1684 1685 split_edges = set([]) # subdivided edges 1686 1687 def triangle(x1, y1, x2, y2, x3, y3, depth=0): 1688 v1, v2, v3 = eval_f(x1, y1), eval_f(x2, y2), eval_f(x3, y3) 1689 1690 if (v1 is v2 is v3 is None) and (depth > max_depth // 2): 1691 # fast finish because the entire region is undefined but 1692 # recurse 'a little' to avoid missing well defined regions 1693 return 1694 elif v1 is None or v2 is None or v3 is None: 1695 # 'triforce' pattern recursion to find the edge of defined region 1696 # 1 1697 # /\ 1698 # 4 /__\ 6 1699 # /\ /\ 1700 # /__\/__\ 1701 # 2 5 3 1702 if depth < max_depth: 1703 x4, y4 = 0.5 * (x1 + x2), 0.5 * (y1 + y2) 1704 x5, y5 = 0.5 * (x2 + x3), 0.5 * (y2 + y3) 1705 x6, y6 = 0.5 * (x1 + x3), 0.5 * (y1 + y3) 1706 split_edges.add( 1707 ((x1, y1), (x2, y2)) 1708 if (x2, y2) > (x1, y1) 1709 else ((x2, y2), (x1, y1)) 1710 ) 1711 split_edges.add( 1712 ((x2, y2), (x3, y3)) 1713 if (x3, y3) > (x2, y2) 1714 else ((x3, y3), (x2, y2)) 1715 ) 1716 split_edges.add( 1717 ((x1, y1), (x3, y3)) 1718 if (x3, y3) > (x1, y1) 1719 else ((x3, y3), (x1, y1)) 1720 ) 1721 triangle(x1, y1, x4, y4, x6, y6, depth + 1) 1722 triangle(x4, y4, x2, y2, x5, y5, depth + 1) 1723 triangle(x6, y6, x5, y5, x3, y3, depth + 1) 1724 triangle(x4, y4, x5, y5, x6, y6, depth + 1) 1725 return 1726 triangles.append(sorted(((x1, y1, v1), (x2, y2, v2), (x3, y3, v3)))) 1727 1728 # linear (grid) sampling 1729 numx = plotpoints[0] * 1.0 1730 numy = plotpoints[1] * 1.0 1731 for xi in range(plotpoints[0]): 1732 for yi in range(plotpoints[1]): 1733 # Decide which way to break the square grid into triangles 1734 # by looking at diagonal lengths. 1735 # 1736 # 3___4 3___4 1737 # |\ | | /| 1738 # | \ | versus | / | 1739 # |__\| |/__| 1740 # 1 2 1 2 1741 # 1742 # Approaching the boundary of the well defined region is 1743 # important too. Use first stategy if 1 or 4 are undefined 1744 # and stategy 2 if either 2 or 3 are undefined. 1745 # 1746 (x1, x2, x3, x4) = ( 1747 xstart + value * (xstop - xstart) 1748 for value in ( 1749 xi / numx, 1750 (xi + 1) / numx, 1751 xi / numx, 1752 (xi + 1) / numx, 1753 ) 1754 ) 1755 (y1, y2, y3, y4) = ( 1756 ystart + value * (ystop - ystart) 1757 for value in ( 1758 yi / numy, 1759 yi / numy, 1760 (yi + 1) / numy, 1761 (yi + 1) / numy, 1762 ) 1763 ) 1764 1765 v1 = eval_f(x1, y1) 1766 v2 = eval_f(x2, y2) 1767 v3 = eval_f(x3, y3) 1768 v4 = eval_f(x4, y4) 1769 1770 if v1 is None or v4 is None: 1771 triangle(x1, y1, x2, y2, x3, y3) 1772 triangle(x4, y4, x3, y3, x2, y2) 1773 elif v2 is None or v3 is None: 1774 triangle(x2, y2, x1, y1, x4, y4) 1775 triangle(x3, y3, x4, y4, x1, y1) 1776 else: 1777 if abs(v3 - v2) > abs(v4 - v1): 1778 triangle(x2, y2, x1, y1, x4, y4) 1779 triangle(x3, y3, x4, y4, x1, y1) 1780 else: 1781 triangle(x1, y1, x2, y2, x3, y3) 1782 triangle(x4, y4, x3, y3, x2, y2) 1783 1784 # adaptive resampling 1785 # TODO: optimise this 1786 # Cos of the maximum angle between successive line segments 1787 ang_thresh = cos(20 * pi / 180) 1788 for depth in range(1, max_depth): 1789 needs_removal = set([]) 1790 lent = len(triangles) # number of initial triangles 1791 for i1 in range(lent): 1792 for i2 in range(lent): 1793 # find all edge pairings 1794 if i1 == i2: 1795 continue 1796 t1 = triangles[i1] 1797 t2 = triangles[i2] 1798 1799 edge_pairing = ( 1800 (t1[0], t1[1]) == (t2[0], t2[1]) 1801 or (t1[0], t1[1]) == (t2[1], t2[2]) 1802 or (t1[0], t1[1]) == (t2[0], t2[2]) 1803 or (t1[1], t1[2]) == (t2[0], t2[1]) 1804 or (t1[1], t1[2]) == (t2[1], t2[2]) 1805 or (t1[1], t1[2]) == (t2[0], t2[2]) 1806 or (t1[0], t1[2]) == (t2[0], t2[1]) 1807 or (t1[0], t1[2]) == (t2[1], t2[2]) 1808 or (t1[0], t1[2]) == (t2[0], t2[2]) 1809 ) 1810 if not edge_pairing: 1811 continue 1812 v1 = [t1[1][i] - t1[0][i] for i in range(3)] 1813 w1 = [t1[2][i] - t1[0][i] for i in range(3)] 1814 v2 = [t2[1][i] - t2[0][i] for i in range(3)] 1815 w2 = [t2[2][i] - t2[0][i] for i in range(3)] 1816 n1 = ( # surface normal for t1 1817 (v1[1] * w1[2]) - (v1[2] * w1[1]), 1818 (v1[2] * w1[0]) - (v1[0] * w1[2]), 1819 (v1[0] * w1[1]) - (v1[1] * w1[0]), 1820 ) 1821 n2 = ( # surface normal for t2 1822 (v2[1] * w2[2]) - (v2[2] * w2[1]), 1823 (v2[2] * w2[0]) - (v2[0] * w2[2]), 1824 (v2[0] * w2[1]) - (v2[1] * w2[0]), 1825 ) 1826 try: 1827 angle = ( 1828 n1[0] * n2[0] + n1[1] * n2[1] + n1[2] * n2[2] 1829 ) / sqrt( 1830 (n1[0] ** 2 + n1[1] ** 2 + n1[2] ** 2) 1831 * (n2[0] ** 2 + n2[1] ** 2 + n2[2] ** 2) 1832 ) 1833 except ZeroDivisionError: 1834 angle = 0.0 1835 if abs(angle) < ang_thresh: 1836 for i, t in ((i1, t1), (i2, t2)): 1837 # subdivide 1838 x1, y1 = t[0][0], t[0][1] 1839 x2, y2 = t[1][0], t[1][1] 1840 x3, y3 = t[2][0], t[2][1] 1841 x4, y4 = 0.5 * (x1 + x2), 0.5 * (y1 + y2) 1842 x5, y5 = 0.5 * (x2 + x3), 0.5 * (y2 + y3) 1843 x6, y6 = 0.5 * (x1 + x3), 0.5 * (y1 + y3) 1844 needs_removal.add(i) 1845 split_edges.add( 1846 ((x1, y1), (x2, y2)) 1847 if (x2, y2) > (x1, y1) 1848 else ((x2, y2), (x1, y1)) 1849 ) 1850 split_edges.add( 1851 ((x2, y2), (x3, y3)) 1852 if (x3, y3) > (x2, y2) 1853 else ((x3, y3), (x2, y2)) 1854 ) 1855 split_edges.add( 1856 ((x1, y1), (x3, y3)) 1857 if (x3, y3) > (x1, y1) 1858 else ((x3, y3), (x1, y1)) 1859 ) 1860 triangle(x1, y1, x4, y4, x6, y6, depth=depth) 1861 triangle(x2, y2, x4, y4, x5, y5, depth=depth) 1862 triangle(x3, y3, x5, y5, x6, y6, depth=depth) 1863 triangle(x4, y4, x5, y5, x6, y6, depth=depth) 1864 # remove subdivided triangles which have been divided 1865 triangles = [ 1866 t for i, t in enumerate(triangles) if i not in needs_removal 1867 ] 1868 1869 # fix up subdivided edges 1870 # 1871 # look at every triangle and see if its edges need updating. 1872 # depending on how many edges require subdivision we proceede with 1873 # one of two subdivision strategies 1874 # 1875 # TODO possible optimisation: don't look at every triangle again 1876 made_changes = True 1877 while made_changes: 1878 made_changes = False 1879 new_triangles = [] 1880 for i, t in enumerate(triangles): 1881 new_points = [] 1882 if ((t[0][0], t[0][1]), (t[1][0], t[1][1])) in split_edges: 1883 new_points.append([0, 1]) 1884 if ((t[1][0], t[1][1]), (t[2][0], t[2][1])) in split_edges: 1885 new_points.append([1, 2]) 1886 if ((t[0][0], t[0][1]), (t[2][0], t[2][1])) in split_edges: 1887 new_points.append([0, 2]) 1888 1889 if len(new_points) == 0: 1890 continue 1891 made_changes = True 1892 # 'triforce' subdivision 1893 # 1 1894 # /\ 1895 # 4 /__\ 6 1896 # /\ /\ 1897 # /__\/__\ 1898 # 2 5 3 1899 # if less than three edges require subdivision bisect them 1900 # anyway but fake their values by averaging 1901 x4 = 0.5 * (t[0][0] + t[1][0]) 1902 y4 = 0.5 * (t[0][1] + t[1][1]) 1903 v4 = stored.get((x4, y4), 0.5 * (t[0][2] + t[1][2])) 1904 1905 x5 = 0.5 * (t[1][0] + t[2][0]) 1906 y5 = 0.5 * (t[1][1] + t[2][1]) 1907 v5 = stored.get((x5, y5), 0.5 * (t[1][2] + t[2][2])) 1908 1909 x6 = 0.5 * (t[0][0] + t[2][0]) 1910 y6 = 0.5 * (t[0][1] + t[2][1]) 1911 v6 = stored.get((x6, y6), 0.5 * (t[0][2] + t[2][2])) 1912 1913 if not (v4 is None or v6 is None): 1914 new_triangles.append(sorted((t[0], (x4, y4, v4), (x6, y6, v6)))) 1915 if not (v4 is None or v5 is None): 1916 new_triangles.append(sorted((t[1], (x4, y4, v4), (x5, y5, v5)))) 1917 if not (v5 is None or v6 is None): 1918 new_triangles.append(sorted((t[2], (x5, y5, v5), (x6, y6, v6)))) 1919 if not (v4 is None or v5 is None or v6 is None): 1920 new_triangles.append( 1921 sorted(((x4, y4, v4), (x5, y5, v5), (x6, y6, v6))) 1922 ) 1923 triangles[i] = None 1924 1925 triangles.extend(new_triangles) 1926 triangles = [t for t in triangles if t is not None] 1927 1928 # add the mesh 1929 mesh_points = [] 1930 if mesh == "System`Full": 1931 for xi in range(plotpoints[0] + 1): 1932 xval = xstart + xi / numx * (xstop - xstart) 1933 mesh_row = [] 1934 for yi in range(plotpoints[1] + 1): 1935 yval = ystart + yi / numy * (ystop - ystart) 1936 z = stored[(xval, yval)] 1937 mesh_row.append((xval, yval, z)) 1938 mesh_points.append(mesh_row) 1939 1940 for yi in range(plotpoints[1] + 1): 1941 yval = ystart + yi / numy * (ystop - ystart) 1942 mesh_col = [] 1943 for xi in range(plotpoints[0] + 1): 1944 xval = xstart + xi / numx * (xstop - xstart) 1945 z = stored[(xval, yval)] 1946 mesh_col.append((xval, yval, z)) 1947 mesh_points.append(mesh_col) 1948 1949 # handle edge subdivisions 1950 made_changes = True 1951 while made_changes: 1952 made_changes = False 1953 for mesh_line in mesh_points: 1954 i = 0 1955 while i < len(mesh_line) - 1: 1956 x1, y1, v1 = mesh_line[i] 1957 x2, y2, v2 = mesh_line[i + 1] 1958 key = ( 1959 ((x1, y1), (x2, y2)) 1960 if (x2, y2) > (x1, y1) 1961 else ((x2, y2), (x1, y1)) 1962 ) 1963 if key in split_edges: 1964 x3 = 0.5 * (x1 + x2) 1965 y3 = 0.5 * (y1 + y2) 1966 v3 = stored[(x3, y3)] 1967 mesh_line.insert(i + 1, (x3, y3, v3)) 1968 made_changes = True 1969 i += 1 1970 i += 1 1971 1972 # handle missing regions 1973 old_meshpoints, mesh_points = mesh_points, [] 1974 for mesh_line in old_meshpoints: 1975 mesh_points.extend( 1976 [ 1977 sorted(g) 1978 for k, g in itertools.groupby( 1979 mesh_line, lambda x: x[2] is None 1980 ) 1981 ] 1982 ) 1983 mesh_points = [ 1984 mesh_line 1985 for mesh_line in mesh_points 1986 if not any(x[2] is None for x in mesh_line) 1987 ] 1988 elif mesh == "System`All": 1989 mesh_points = set([]) 1990 for t in triangles: 1991 mesh_points.add((t[0], t[1]) if t[1] > t[0] else (t[1], t[0])) 1992 mesh_points.add((t[1], t[2]) if t[2] > t[1] else (t[2], t[1])) 1993 mesh_points.add((t[0], t[2]) if t[2] > t[0] else (t[2], t[0])) 1994 mesh_points = list(mesh_points) 1995 1996 # find the max and min height 1997 v_min = v_max = None 1998 for t in triangles: 1999 for tx, ty, v in t: 2000 if v_min is None or v < v_min: 2001 v_min = v 2002 if v_max is None or v > v_max: 2003 v_max = v 2004 graphics.extend( 2005 self.construct_graphics( 2006 triangles, mesh_points, v_min, v_max, options, evaluation 2007 ) 2008 ) 2009 return self.final_graphics(graphics, options) 2010 2011 2012class Plot(_Plot): 2013 """ 2014 <dl> 2015 <dt>'Plot[$f$, {$x$, $xmin$, $xmax$}]' 2016 <dd>plots $f$ with $x$ ranging from $xmin$ to $xmax$. 2017 2018 <dt>'Plot[{$f1$, $f2$, ...}, {$x$, $xmin$, $xmax$}]' 2019 <dd>plots several functions $f1$, $f2$, ... 2020 2021 </dl> 2022 2023 >> Plot[{Sin[x], Cos[x], x / 3}, {x, -Pi, Pi}] 2024 = -Graphics- 2025 2026 >> Plot[Sin[x], {x, 0, 4 Pi}, PlotRange->{{0, 4 Pi}, {0, 1.5}}] 2027 = -Graphics- 2028 2029 >> Plot[Tan[x], {x, -6, 6}, Mesh->Full] 2030 = -Graphics- 2031 2032 >> Plot[x^2, {x, -1, 1}, MaxRecursion->5, Mesh->All] 2033 = -Graphics- 2034 2035 >> Plot[Log[x], {x, 0, 5}, MaxRecursion->0] 2036 = -Graphics- 2037 2038 >> Plot[Tan[x], {x, 0, 6}, Mesh->All, PlotRange->{{-1, 5}, {0, 15}}, MaxRecursion->10] 2039 = -Graphics- 2040 2041 A constant function: 2042 >> Plot[3, {x, 0, 1}] 2043 = -Graphics- 2044 2045 #> Plot[1 / x, {x, -1, 1}] 2046 = -Graphics- 2047 #> Plot[x, {y, 0, 2}] 2048 = -Graphics- 2049 2050 #> Plot[{f[x],-49x/12+433/108},{x,-6,6}, PlotRange->{-10,10}, AspectRatio->{1}] 2051 = -Graphics- 2052 2053 #> Plot[Sin[t], {t, 0, 2 Pi}, PlotPoints -> 1] 2054 : Value of option PlotPoints -> 1 is not an integer >= 2. 2055 = Plot[Sin[t], {t, 0, 2 Pi}, PlotPoints -> 1] 2056 2057 #> Plot[x*y, {x, -1, 1}] 2058 = -Graphics- 2059 """ 2060 2061 def get_functions_param(self, functions): 2062 if functions.has_form("List", None): 2063 functions = functions.leaves 2064 else: 2065 functions = [functions] 2066 return functions 2067 2068 def get_plotrange(self, plotrange, start, stop): 2069 x_range = y_range = None 2070 if isinstance(plotrange, numbers.Real): 2071 plotrange = ["System`Full", [-plotrange, plotrange]] 2072 if plotrange == "System`Automatic": 2073 plotrange = ["System`Full", "System`Automatic"] 2074 elif plotrange == "System`All": 2075 plotrange = ["System`All", "System`All"] 2076 if isinstance(plotrange, list) and len(plotrange) == 2: 2077 if isinstance(plotrange[0], numbers.Real) and isinstance( # noqa 2078 plotrange[1], numbers.Real 2079 ): 2080 x_range, y_range = "System`Full", plotrange 2081 else: 2082 x_range, y_range = plotrange 2083 if x_range == "System`Full": 2084 x_range = [start, stop] 2085 return x_range, y_range 2086 2087 def eval_f(self, f, x_value): 2088 value = f(x_value) 2089 if value is not None: 2090 return (x_value, value) 2091 2092 2093class ParametricPlot(_Plot): 2094 """ 2095 <dl> 2096 <dt>'ParametricPlot[{$f_x$, $f_y$}, {$u$, $umin$, $umax$}]' 2097 <dd>plots a parametric function $f$ with the parameter $u$ ranging from $umin$ to $umax$. 2098 <dt>'ParametricPlot[{{$f_x$, $f_y$}, {$g_x$, $g_y$}, ...}, {$u$, $umin$, $umax$}]' 2099 <dd>plots several parametric functions $f$, $g$, ... 2100 <dt>'ParametricPlot[{$f_x$, $f_y$}, {$u$, $umin$, $umax$}, {$v$, $vmin$, $vmax$}]' 2101 <dd>plots a parametric area. 2102 <dt>'ParametricPlot[{{$f_x$, $f_y$}, {$g_x$, $g_y$}, ...}, {$u$, $umin$, $umax$}, {$v$, $vmin$, $vmax$}]' 2103 <dd>plots several parametric areas. 2104 </dl> 2105 2106 >> ParametricPlot[{Sin[u], Cos[3 u]}, {u, 0, 2 Pi}] 2107 = -Graphics- 2108 2109 >> ParametricPlot[{Cos[u] / u, Sin[u] / u}, {u, 0, 50}, PlotRange->0.5] 2110 = -Graphics- 2111 2112 >> ParametricPlot[{{Sin[u], Cos[u]},{0.6 Sin[u], 0.6 Cos[u]}, {0.2 Sin[u], 0.2 Cos[u]}}, {u, 0, 2 Pi}, PlotRange->1, AspectRatio->1] 2113 = -Graphics- 2114 """ 2115 2116 expect_list = True 2117 2118 def get_functions_param(self, functions): 2119 if functions.has_form("List", 2) and not ( 2120 functions.leaves[0].has_form("List", None) 2121 or functions.leaves[1].has_form("List", None) 2122 ): 2123 # One function given 2124 functions = [functions] 2125 else: 2126 # Multiple Functions 2127 functions = functions.leaves 2128 return functions 2129 2130 def get_plotrange(self, plotrange, start, stop): 2131 x_range = y_range = None 2132 if isinstance(plotrange, numbers.Real): 2133 plotrange = [[-plotrange, plotrange], [-plotrange, plotrange]] 2134 if plotrange == "System`Automatic": 2135 plotrange = ["System`Automatic", "System`Automatic"] 2136 elif plotrange == "System`All": 2137 plotrange = ["System`All", "System`All"] 2138 if isinstance(plotrange, list) and len(plotrange) == 2: 2139 if isinstance(plotrange[0], numbers.Real) and isinstance( # noqa 2140 plotrange[1], numbers.Real 2141 ): 2142 x_range = [-plotrange[0], plotrange[1]] 2143 y_range = [-plotrange[1], plotrange[1]] 2144 else: 2145 x_range, y_range = plotrange 2146 return x_range, y_range 2147 2148 def eval_f(self, f, x_value): 2149 value = f(x_value) 2150 if value is not None and len(value) == 2: 2151 return value 2152 2153 2154class PolarPlot(_Plot): 2155 """ 2156 <dl> 2157 <dt>'PolarPlot[$r$, {$t$, $tmin$, $tmax$}]' 2158 <dd>creates a polar plot of $r$ with angle $t$ ranging from 2159 $tmin$ to $tmax$. 2160 </dl> 2161 2162 >> PolarPlot[Cos[5t], {t, 0, Pi}] 2163 = -Graphics- 2164 2165 >> PolarPlot[{1, 1 + Sin[20 t] / 5}, {t, 0, 2 Pi}] 2166 = -Graphics- 2167 """ 2168 2169 options = _Plot.options.copy() 2170 options.update( 2171 {"AspectRatio": "1",} 2172 ) 2173 2174 def get_functions_param(self, functions): 2175 if functions.has_form("List", None): 2176 functions = functions.leaves 2177 else: 2178 functions = [functions] 2179 return functions 2180 2181 def get_plotrange(self, plotrange, start, stop): 2182 x_range = y_range = None 2183 if isinstance(plotrange, numbers.Real): 2184 plotrange = [[-plotrange, plotrange], [-plotrange, plotrange]] 2185 if plotrange == "System`Automatic": 2186 plotrange = ["System`Automatic", "System`Automatic"] 2187 elif plotrange == "System`All": 2188 plotrange = ["System`All", "System`All"] 2189 if isinstance(plotrange, list) and len(plotrange) == 2: 2190 if isinstance(plotrange[0], numbers.Real) and isinstance( # noqa 2191 plotrange[1], numbers.Real 2192 ): 2193 x_range = [-plotrange[0], plotrange[1]] 2194 y_range = [-plotrange[1], plotrange[1]] 2195 else: 2196 x_range, y_range = plotrange 2197 return x_range, y_range 2198 2199 def eval_f(self, f, x_value): 2200 value = f(x_value) 2201 if value is not None: 2202 return (value * cos(x_value), value * sin(x_value)) 2203 2204 2205class ListPlot(_ListPlot): 2206 """ 2207 <dl> 2208 <dt>'ListPlot[{$y_1$, $y_2$, ...}]' 2209 <dd>plots a list of y-values, assuming integer x-values 1, 2, 3, ... 2210 <dt>'ListPlot[{{$x_1$, $y_1$}, {$x_2$, $y_2$}, ...}]' 2211 <dd>plots a list of $x$, $y$ pairs. 2212 <dt>'ListPlot[{$list_1$, $list_2$, ...}]' 2213 <dd>plots several lists of points. 2214 </dl> 2215 2216 ListPlot accepts a superset of the Graphics options. 2217 2218 >> ListPlot[Table[n ^ 2, {n, 10}]] 2219 = -Graphics- 2220 """ 2221 2222 attributes = ("HoldAll",) 2223 2224 options = Graphics.options.copy() 2225 options.update( 2226 { 2227 "Axes": "True", 2228 "AspectRatio": "1 / GoldenRatio", 2229 "Mesh": "None", 2230 "PlotRange": "Automatic", 2231 "PlotPoints": "None", 2232 "Filling": "None", 2233 "Joined": "False", 2234 } 2235 ) 2236 2237 2238class ListLinePlot(_ListPlot): 2239 """ 2240 <dl> 2241 <dt>'ListLinePlot[{$y_1$, $y_2$, ...}]' 2242 <dd>plots a line through a list of $y$-values, assuming integer $x$-values 1, 2, 3, ... 2243 2244 <dt>'ListLinePlot[{{$x_1$, $y_1$}, {$x_2$, $y_2$}, ...}]' 2245 <dd>plots a line through a list of $x$, $y$ pairs. 2246 2247 <dt>'ListLinePlot[{$list_1$, $list_2$, ...}]' 2248 <dd>plots several lines. 2249 </dl> 2250 2251 ListPlot accepts a superset of the Graphics options. 2252 2253 >> ListLinePlot[Table[{n, n ^ 0.5}, {n, 10}]] 2254 = -Graphics- 2255 2256 >> ListLinePlot[{{-2, -1}, {-1, -1}}] 2257 = -Graphics- 2258 """ 2259 2260 attributes = ("HoldAll",) 2261 2262 options = Graphics.options.copy() 2263 options.update( 2264 { 2265 "Axes": "True", 2266 "AspectRatio": "1 / GoldenRatio", 2267 "Mesh": "None", 2268 "PlotRange": "Automatic", 2269 "PlotPoints": "None", 2270 "Filling": "None", 2271 "Joined": "True", 2272 } 2273 ) 2274 2275 2276class Plot3D(_Plot3D): 2277 """ 2278 <dl> 2279 <dt>'Plot3D[$f$, {$x$, $xmin$, $xmax$}, {$y$, $ymin$, $ymax$}]' 2280 <dd>creates a three-dimensional plot of $f$ with $x$ ranging from $xmin$ to $xmax$ and $y$ ranging from $ymin$ to $ymax$. 2281 2282 </dl> 2283 2284 Plot3D has the same options as Graphics3D, in particular: 2285 <ul> 2286 <li>Mesh 2287 <li>PlotPoints 2288 <li>MaxRecursion 2289 </ul> 2290 2291 2292 >> Plot3D[x ^ 2 + 1 / y, {x, -1, 1}, {y, 1, 4}] 2293 = -Graphics3D- 2294 2295 >> Plot3D[Sin[y + Sin[3 x]], {x, -2, 2}, {y, -2, 2}, PlotPoints->20] 2296 = -Graphics3D- 2297 2298 >> Plot3D[x / (x ^ 2 + y ^ 2 + 1), {x, -2, 2}, {y, -2, 2}, Mesh->None] 2299 = -Graphics3D- 2300 2301 >> Plot3D[Sin[x y] /(x y), {x, -3, 3}, {y, -3, 3}, Mesh->All] 2302 = -Graphics3D- 2303 2304 >> Plot3D[Log[x + y^2], {x, -1, 1}, {y, -1, 1}] 2305 = -Graphics3D- 2306 2307 #> Plot3D[z, {x, 1, 20}, {y, 1, 10}] 2308 = -Graphics3D- 2309 2310 ## MaxRecursion Option 2311 #> Plot3D[0, {x, -2, 2}, {y, -2, 2}, MaxRecursion -> 0] 2312 = -Graphics3D- 2313 #> Plot3D[0, {x, -2, 2}, {y, -2, 2}, MaxRecursion -> 15] 2314 = -Graphics3D- 2315 #> Plot3D[0, {x, -2, 2}, {y, -2, 2}, MaxRecursion -> 16] 2316 : MaxRecursion must be a non-negative integer; the recursion value is limited to 15. Using MaxRecursion -> 15. 2317 = -Graphics3D- 2318 #> Plot3D[0, {x, -2, 2}, {y, -2, 2}, MaxRecursion -> -1] 2319 : MaxRecursion must be a non-negative integer; the recursion value is limited to 15. Using MaxRecursion -> 0. 2320 = -Graphics3D- 2321 #> Plot3D[0, {x, -2, 2}, {y, -2, 2}, MaxRecursion -> a] 2322 : MaxRecursion must be a non-negative integer; the recursion value is limited to 15. Using MaxRecursion -> 0. 2323 = -Graphics3D- 2324 #> Plot3D[0, {x, -2, 2}, {y, -2, 2}, MaxRecursion -> Infinity] 2325 : MaxRecursion must be a non-negative integer; the recursion value is limited to 15. Using MaxRecursion -> 15. 2326 = -Graphics3D- 2327 2328 #> Plot3D[x ^ 2 + 1 / y, {x, -1, 1}, {y, 1, z}] 2329 : Limiting value z in {y, 1, z} is not a machine-size real number. 2330 = Plot3D[x ^ 2 + 1 / y, {x, -1, 1}, {y, 1, z}] 2331 """ 2332 2333 # FIXME: This test passes but the result is 511 lines long ! 2334 """ 2335 #> Plot3D[x + 2y, {x, -2, 2}, {y, -2, 2}] // TeXForm 2336 """ 2337 2338 attributes = ("HoldAll",) 2339 2340 options = Graphics.options.copy() 2341 options.update( 2342 { 2343 "Axes": "True", 2344 "AspectRatio": "1", 2345 "Mesh": "Full", 2346 "PlotPoints": "None", 2347 "BoxRatios": "{1, 1, 0.4}", 2348 "MaxRecursion": "2", 2349 } 2350 ) 2351 2352 def get_functions_param(self, functions): 2353 if functions.has_form("List", None): 2354 return functions.leaves 2355 else: 2356 return [functions] 2357 2358 def construct_graphics( 2359 self, triangles, mesh_points, v_min, v_max, options, evaluation 2360 ): 2361 graphics = [] 2362 for p1, p2, p3 in triangles: 2363 graphics.append( 2364 Expression( 2365 "Polygon", 2366 Expression( 2367 "List", 2368 Expression(SymbolList, *p1), 2369 Expression(SymbolList, *p2), 2370 Expression(SymbolList, *p3), 2371 ), 2372 ) 2373 ) 2374 # Add the Grid 2375 for xi in range(len(mesh_points)): 2376 line = [] 2377 for yi in range(len(mesh_points[xi])): 2378 line.append( 2379 Expression( 2380 "List", 2381 mesh_points[xi][yi][0], 2382 mesh_points[xi][yi][1], 2383 mesh_points[xi][yi][2], 2384 ) 2385 ) 2386 graphics.append(Expression("Line", Expression(SymbolList, *line))) 2387 return graphics 2388 2389 def final_graphics(self, graphics, options): 2390 return Expression( 2391 "Graphics3D", 2392 Expression(SymbolList, *graphics), 2393 *options_to_rules(options, Graphics3D.options) 2394 ) 2395 2396 2397class DensityPlot(_Plot3D): 2398 """ 2399 <dl> 2400 <dt>'DensityPlot[$f$, {$x$, $xmin$, $xmax$}, {$y$, $ymin$, $ymax$}]' 2401 <dd>plots a density plot of $f$ with $x$ ranging from $xmin$ to $xmax$ and $y$ ranging from $ymin$ to $ymax$. 2402 </dl> 2403 2404 >> DensityPlot[x ^ 2 + 1 / y, {x, -1, 1}, {y, 1, 4}] 2405 = -Graphics- 2406 2407 >> DensityPlot[1 / x, {x, 0, 1}, {y, 0, 1}] 2408 = -Graphics- 2409 2410 >> DensityPlot[Sqrt[x * y], {x, -1, 1}, {y, -1, 1}] 2411 = -Graphics- 2412 2413 >> DensityPlot[1/(x^2 + y^2 + 1), {x, -1, 1}, {y, -2,2}, Mesh->Full] 2414 = -Graphics- 2415 2416 >> DensityPlot[x^2 y, {x, -1, 1}, {y, -1, 1}, Mesh->All] 2417 = -Graphics- 2418 """ 2419 2420 attributes = ("HoldAll",) 2421 2422 options = Graphics.options.copy() 2423 options.update( 2424 { 2425 "Axes": "False", 2426 "AspectRatio": "1", 2427 "Mesh": "None", 2428 "Frame": "True", 2429 "ColorFunction": "Automatic", 2430 "ColorFunctionScaling": "True", 2431 "PlotPoints": "None", 2432 "MaxRecursion": "0", 2433 # 'MaxRecursion': '2', # FIXME causes bugs in svg output see #303 2434 } 2435 ) 2436 2437 def get_functions_param(self, functions): 2438 return [functions] 2439 2440 def construct_graphics( 2441 self, triangles, mesh_points, v_min, v_max, options, evaluation 2442 ): 2443 color_function = self.get_option(options, "ColorFunction", evaluation, pop=True) 2444 color_function_scaling = self.get_option( 2445 options, "ColorFunctionScaling", evaluation, pop=True 2446 ) 2447 2448 color_function_min = color_function_max = None 2449 if color_function.get_name() == "System`Automatic": 2450 color_function = String("LakeColors") 2451 if color_function.get_string_value(): 2452 func = Expression("ColorData", color_function.get_string_value()).evaluate( 2453 evaluation 2454 ) 2455 if func.has_form("ColorDataFunction", 4): 2456 color_function_min = func.leaves[2].leaves[0].round_to_float() 2457 color_function_max = func.leaves[2].leaves[1].round_to_float() 2458 color_function = Expression( 2459 "Function", Expression(func.leaves[3], Expression("Slot", 1)) 2460 ) 2461 else: 2462 evaluation.message("DensityPlot", "color", func) 2463 return 2464 if color_function.has_form("ColorDataFunction", 4): 2465 color_function_min = color_function.leaves[2].leaves[0].round_to_float() 2466 color_function_max = color_function.leaves[2].leaves[1].round_to_float() 2467 2468 color_function_scaling = color_function_scaling.is_true() 2469 v_range = v_max - v_min 2470 2471 if v_range == 0: 2472 v_range = 1 2473 2474 if color_function.has_form("ColorDataFunction", 4): 2475 color_func = color_function.leaves[3] 2476 else: 2477 color_func = color_function 2478 if ( 2479 color_function_scaling 2480 and color_function_min is not None # noqa 2481 and color_function_max is not None 2482 ): 2483 color_function_range = color_function_max - color_function_min 2484 2485 colors = {} 2486 2487 def eval_color(x, y, v): 2488 v_scaled = (v - v_min) / v_range 2489 if ( 2490 color_function_scaling 2491 and color_function_min is not None # noqa 2492 and color_function_max is not None 2493 ): 2494 v_color_scaled = color_function_min + v_scaled * color_function_range 2495 else: 2496 v_color_scaled = v 2497 2498 # Calculate and store 100 different shades max. 2499 v_lookup = int(v_scaled * 100 + 0.5) 2500 2501 value = colors.get(v_lookup) 2502 if value is None: 2503 value = Expression(color_func, Real(v_color_scaled)) 2504 value = value.evaluate(evaluation) 2505 colors[v_lookup] = value 2506 return value 2507 2508 points = [] 2509 vertex_colors = [] 2510 graphics = [] 2511 for p in triangles: 2512 points.append( 2513 Expression(SymbolList, *(Expression(SymbolList, *x[:2]) for x in p)) 2514 ) 2515 vertex_colors.append(Expression(SymbolList, *(eval_color(*x) for x in p))) 2516 2517 graphics.append( 2518 Expression( 2519 "Polygon", 2520 Expression(SymbolList, *points), 2521 Expression( 2522 "Rule", 2523 Symbol("VertexColors"), 2524 Expression(SymbolList, *vertex_colors), 2525 ), 2526 ) 2527 ) 2528 2529 # add mesh 2530 for xi in range(len(mesh_points)): 2531 line = [] 2532 for yi in range(len(mesh_points[xi])): 2533 line.append( 2534 Expression( 2535 SymbolList, mesh_points[xi][yi][0], mesh_points[xi][yi][1] 2536 ) 2537 ) 2538 graphics.append(Expression("Line", Expression(SymbolList, *line))) 2539 2540 return graphics 2541 2542 def final_graphics(self, graphics, options): 2543 return Expression( 2544 "Graphics", 2545 Expression(SymbolList, *graphics), 2546 *options_to_rules(options, Graphics.options) 2547 ) 2548