Plotting

How can we plot mathematical objects in sagemath?

A directed graph consists of vertices and edges. Edges are directed, thus have a source and destination. We specify an edge as some id (e.g. 1) associated with a list of ids it points to (e.g. [2, 3]):

In [1]:
d = {1: [2, 3], 2: [4], 3: [5, 6], 4: [7], 5: [7], 6: [6, 7], 7: [8]}
g = DiGraph(d)
g.plot()
Out[1]:

One typical problem is shown here: we plot a function on a given domain (here: from ‘-5’ to ‘5’), but the function is only defined on a smaller domain (namely ‘0’ to ‘5’):

In [2]:
plot(sqrt(x), x, -5, 5)

verbose 0 (3763: plot.py, generate_plot_points) WARNING: When plotting, failed to evaluate function at 100 points.
verbose 0 (3763: plot.py, generate_plot_points) Last error message: 'math domain error'
Out[2]:

We can provide keyword-argument ‘detect_poles’, but this will often fail (as pole detection is mathematically difficult). However, I present one function which shows the functionality nicely:

In [3]:
plot(tan(2*x), x, -pi, pi, ymin=-5, ymax=5, detect_poles='show')
Out[3]:

The actual solution is to define a piecewise function (mathematically this extends the domain of the partial function). Here, we slice the domain into to intervals [-5, 0] and (0, 5). Recognize that lists are used to denote closed intervals (matches mathematical notation) and tuples are used for open intervals (matches mathematical notation). ‘0’ denotes the constant zero function:

In [4]:
f = piecewise([
    ([-5, 0], 0),       # f(x) = 0        on [-5, 0]
    ((0, 5), sqrt(x))   # f(x) = sqrt(x)  on (0, 5)
])

plot(f)
Out[4]:

How about mixed intervals? Use the RealSet methods to retrieve subsets of the real line:

In [5]:
f = piecewise([
    ([-5, 0], 0),                         # f(x) = 0        on [-5, 0]
    (RealSet.open_closed(0,5), sqrt(x))   # f(x) = sqrt(x)  on (0, 5]
])

plot(f)
Out[5]:

3D plots also look cool. Here a function takes x- and y-values and returns the corresponding z-value:

In [6]:
f(x, y) = sin(x) * cos(y)
plot3d(f, (x, -pi, pi), (y, -pi, pi))
Out[6]:

Parametric plots are also possible. In case of parametric plots, the input argument(s) is not a dimension itself, but iterates over a domain. The given function returns a tuple specifying all values in all dimensions. Here Neil's parabola is used as an example:

In [7]:
var('x')
parametric_plot((x^2, x^3), (x, -1, 1))
Out[7]:

One super-cool plot from the documentation are these interwoven tori:

In [8]:
# source: https://doc.sagemath.org/html/en/reference/plot3d/sage/plot/plot3d/parametric_plot3d.html#sage.plot.plot3d.parametric_plot3d.parametric_plot3d
u, v = var('u,v')
f1 = (4 + (3 + cos(v)) * sin(u), 4 + (3 + cos(v)) * cos(u), 4 + sin(v))
f2 = (8 + (3 + cos(v)) * cos(u), 3 + sin(v), 4 + (3 + cos(v)) * sin(u))
p1 = parametric_plot3d(f1, (u, 0, 2 * pi), (v, 0, 2 * pi), texture="red")
p2 = parametric_plot3d(f2, (u, 0, 2 * pi), (v, 0, 2 * pi), texture="blue")
p1 + p2
Out[8]:

This polar(-coordinate) plot from the documentation also looks super-nice:

In [9]:
# http://doc.sagemath.org/html/en/reference/plotting/sage/plot/plot.html#sage.plot.plot.polar_plot
polar_plot([(1.2 + k * 0.2) * log(x) for k in range(6)], 1, 3 * pi, fill={0: [1], 2: [3], 4: [5]})
Out[9]:

Finally, you can also make bar charts, but recognize the issues:

  • the legend shows too narrow color boxes
  • It is not possible to annotate all bars along the x-axis

Personally, I concluded that this is one area where sagemath is impractical in terms of plotting.

In [10]:
# data source: deleted 'Grazer Linuxtage' DE-Wikipedia article
Unknown = 0
lectures = [17, 31, 23, 28, 29, 29, 29, 29, 35, 34, 27, 33, 40, 42, 35, 37, 45]
workshops = [Unknown, 3, 6, 5, 4, Unknown, Unknown, 1, 1, 5, 3, 12, 13, 12, 8, 12, 17]

fig = plot([], title="Number of accepted GLT submissions", axes_labels=['years', '#'])
fig += bar_chart(lectures, color='r', legend_label='lectures')
fig += bar_chart(workshops, color='b', legend_label='workshops')
fig.show(show_legend=True, legend_loc="upper center")

Finally, we can also draw simple line segments:

In [11]:
g = Graphics()
g += line([(0, 0), (0, 1), (1, 1), (1, 2), (2, 2)])
g.show()

Giving rise to the following Dragon curve implementation:

In [12]:
Left = 1
Right = 0

iterations = 5
new_direction = Left
coords = [(0, 0), (1, 0)]
colors = '#9400D3;#4B0082;#0000FF;#00FF00;#DDCC00;#FF7F00;#FF0000'.split(';')

def reverse_number(num):
    """Reverse bits of a number; e.g. LSB becomes MSB"""
    num = int(num)
    b = num.bit_length()
    result = int(0)
    for i in range(b):
        result |= ((num >> i) & 1) << (b - i - 1)
    return int(result)

assert reverse_number(0b111001) == 0b100111

def iterate(turns):
    """Iterate `turns`. Let `abc` be the digits of a binary representation of turns less than 8.
    Then the result is `abcXcÌ…bÌ…aÌ…` where `X` is 1 if new_direction is Left or 0 if new_direction is Right
    and where `cÌ…bÌ…aÌ…` represents `cba` with each bit inverted.
    """
    turns = int(turns)
    b = turns.bit_length()
    return (turns << (b + 1)) | (new_direction << b) | (~reverse_number(turns) & (2**b - 1))

assert new_direction == Left and Left == 1
assert iterate(0b1100) == 0b110011100

def next_point(a, b, turn_direction=Left):
    """Given a 2-dimensional point in a Cartesian coordinate system as tuple,
    find the next point if we continue at point `b` in direction `turn_direction` with length `b - a`.

    :param a:               last point
    :param b:               previous to last point
    :param turn_direction:  direction to turn
    """
    factor = -1 if (turn_direction == Left) else 1
    return (b[0] + factor * (b[1] - a[1]),
            b[1] - factor * (b[0] - a[0]))

assert next_point((3, 3), (3, 2), True) == (4, 2)
assert next_point((3, 3), (3, 2), False) == (2, 2)

def plot_lines(fig, coords, color):
    """Given a list of (x,y) coordinates, generate an image """
    xmin = min(x for (x, y) in coords)
    xmax = max(x for (x, y) in coords)
    ymin = min(y for (x, y) in coords)
    ymax = max(y for (x, y) in coords)
    
    assert xmin <= xmax
    assert ymin <= ymax

    fig += line(coords, rgbcolor=color)

    # proper parameters for the plot
    mm = fig.get_minmax_data()
    fig.xmin(min(xmin, mm['xmin']))
    fig.xmax(max(xmax, mm['xmax']))
    fig.ymin(min(ymin, mm['ymin']))
    fig.ymax(max(ymax, mm['ymax']))
    fig.set_aspect_ratio(1)
    fig.axes(False)
    
    return fig

def main(initial_pos=None, iterations=5):
    """Main routine"""
    # parse/prepare arguments
    if initial_pos is None:
        coords = [(0, 0), (1, 0)]
    else:
        coords = initial_pos

    # compute turns we need to do
    fig = Graphics()
    for i in range(iterations):
        turns = 0
        for _ in range(i):
            turns = iterate(turns)
        #print('turns: {}'.format(bin(turns)[2:].replace(str(Right), 'R').replace(str(Left), 'L')))

        # compute next coordinate
        coords = coords[0:2]
        for j in range(2**i - 1):
            direction = bool((turns >> (2**i - j - 2)) & 1)
            #print('Left' if direction == Left else 'Right')
            next_pt = next_point(coords[-2], coords[-1], direction)
            coords.append(next_pt)
        #print('coordinates: {}'.format(coords))

        # plot
        fig = plot_lines(fig, coords[len(coords) // 2:], colors[i % len(colors)])

    show(fig)
    fig.save('dragon_curve.svg')

main(iterations=13)
In [ ]: