Revision as of 22:28, 17 February 2021

Crash Course in Scientific Computing

V. Plotting

Last lecture we have left the Conway's game of Life cellular automaton computing in arrays which were not very easy to visualize. Actually, we could have made a better job at displaying them for instance by using better characters to represent live and dead cells:

include("life.jl")
cleanspace(12);
putglider(6,6);
replace(space, 1=>*, 0=>-)

12×12 Array{Any,2}:
 -  -  -  -  -  -  -  -  -  -  -  -
 -  -  -  -  -  -  -  -  -  -  -  -
 -  -  -  -  -  -  -  -  -  -  -  -
 -  -  -  -  *  -  -  -  -  -  -  -
 -  -  -  -  -  *  -  -  -  -  -  -
 -  -  -  *  *  *  -  -  -  -  -  -
 -  -  -  -  -  -  -  -  -  -  -  -
 -  -  -  -  -  -  -  -  -  -  -  -
 -  -  -  -  -  -  -  -  -  -  -  -
 -  -  -  -  -  -  -  -  -  -  -  -
 -  -  -  -  -  -  -  -  -  -  -  -
 -  -  -  -  -  -  -  -  -  -  -  -

But a computer can do much better than that, by using computer graphics. This is an ample and complex problem, of which one can only hope to scratch the surface. That'll be providing plenty of material already.

For such a lecture deeply rooted on the graphical output, we will be using the notebook feature. This will allows us to keep all our graphs since the REPL would update them in a separate window. Jupyter (the notebook platform) has several shortcuts which it is handy to keep in mind when editing:

• One can toggle between edit (green) and command (blue) mode with esc and enter respectively.
• Use a and b to insert a cell above and below the current position.
• Use d d to delete the current cell
• Use z to undo
• Use m to transform a cell into markdown, which is the non-code part of the notebook that documents it. Further commands of the markdown are # for a section, ## for a subsection, etc. One can also use LaTeX there. We shall see some examples.

Plotting is not native in Julia, which you might regard as a pity, since it's such a fundamental aspect of computer programming. But this is actually a feature, as it allows us to choose the plotting environment that we like. We are not limited to built-in plotting capacities. We will be using Plots:

using Plots

This is a plotting metapackage which brings many different plotting packages under a single API, including Plotly, which is a powerful library used predominantly by Python.

So let's come back to our cellular automaton. Let's use a bigger space (ASCII limitations won't be a problem now as we'll be able to get down to the pixel level):

szspace=100;
cleanspace(szspace);
putglider(szspace÷2,szspace÷2,) # put a glider in the middle

We can then use heatmap() to make a bitmap rendering of the array. In the notebook environment, this gives us the following:

We now can better see the pattern which was so difficult to discern when scrambled in a table of 0 and 1. The default plotting is still far from satisfactory since colors are not particularly pretty, we don't need the barscale and the aspect ratio is favouring one dimension, with no reason. All this can be remedied with options:

heatmap(space,c=cgrad([:white, :black]), axis=nothing, legend=:none, aspect_ratio=1, yflip=true)

This will be best encapsulated in a function:

function plotspace(space)
    heatmap(transpose(space), c=cgrad([:white, :black]), axis=nothing, legend=:none, aspectratio=1, yflip=true)
end

Note that we transpose the space and flip the $y$ axis, so as to have the display read with the same conventions as the matrix (i.e., with (1,1) on the top left and (szspace, szspace) at the bottom right).

We can now let our Conway universe evolve: .

nextgen();plotspace(space)

which reveals the surprise first discovered by Richard Guy, that this little pattern replicates itself at the same time as it translates in the space. Things can move! This was actually discovered by following the evolution of the so-called r-pentomino, which was interesting Conway and his group for its capacity to live a very long time (this was one appeal of the system, to study its capacity to evolve in an irregular way for very long times). The evolution of the R-pentomino is worth visualizing. Like this, we can replicate this epic discovery of computer science:

cleanspace(szspace);
putrpentomino(szspace÷2,szspace÷2)

rpento = @animate for i ∈ 1:200
           nextgen();
           plotspace(space)
       end
gif(rpento, fps=5)

The irresistible movement of the glider is now clear on such an animation. This is a major result of computer science. The rules of the automaton are extremely simple: they bear no notion of objects beside that of a single-cell, and even less of notions such as mutation, movement, interactions, etc., and all this emerges all by itself from the rules. This led some people to postulate that the laws of Physics are similarly emerging from very basic rules. We will come back to that. For now, our emphasis being on graphical features, we will

gosper = @animate for i ∈ 1:1000
           nextgen();
           heatmap(space,c=cgrad([:white, :black]), axis=nothing, legend=:none, aspect_ratio=1, yflip=true)
       end
gif(gosper, fps=25)

From this point onward, we can literally seat and see: there is much that can be investigated, for instance, we could aim two guns at each other and see how they interact. This is one example out of many:

szspace=120;
cleanspace(szspace);
putGosperglidergun(2,2);
offsetx=29;
offsety=35;
space[szspace-offsetx:szspace-offsetx+9-1, szspace-offsety:szspace-offsety+36-1]=reverse(reverse(Gosperglidergun(),dims=1),dims=2)
gosper = @animate for i ∈ 1:500
           nextgen();
           plotspace(space)
       end
gif(gosper, fps=50)

We let you explore further, if you're puzzled by this simple universe, and will return to the study of ours.

Plotting bitmaps is of course not the most common application we want from a computer. Maybe the most wanted overall is to plot raw data, which typically comes as a list of points or a list of pair of points (to be interpreted in the $(x, y)$ or $y_1, y_2$ sense, making one or two curves, respectively). Let us look at examples. For instance, let us track a Collatz trajectory, for which we'll redefine slightly our collatz function so that it gives us a list of the elements:

function collatz(n)
    global lcol=[n];
    while n!=1
        if n%2==0
            n÷=2
        else
            n=3n+1
        end
        push!(lcol, n)
    end
    lcol
end

We have used push! to append the element to the list, so that we don't need to predefine arrays, given that we don't know their size ahead of time. There is also an append command which allows to append more than one element. Now we can see how the series evolves:

plot(collatz(3), lw=2, xlabel="n", ylabel="collatz series", label="3")

We can use plot! to add to the previous result. Note that in Julia, ! refers to a function that alters its argument (we've seen that with push! already):

plot!(collatz(5), lw=2, xlabel="n", ylabel="collatz series", label="5")

Like this it's easy to combine plots and computer programming:

p = plot()
for k ∈ 2:10
    global p=plot!(collatz(k), lw=2, xlabel="n", ylabel="collatz series", label=k)
end
p

There is a lot of redundancy in this graph as various series collapse onto the same cycle. That makes for an interesting exercise for a better visualization of the Collatz series.

Maybe the trajectories of each series is not the most important, enlightening or interesting, baring exceptions:

We would more interesting in, say, the lengths. Let's change our collatz function again:

function collatz(n)
    global counter=1;
    while n!=1
        if n%2==0
            n÷=2
        else
            n=3n+1
        end
        counter+=1;
    end
    counter
end

In which case it is trivial to plot the lengths:

plot([collatz(i) for i ∈ 1:100], xlabel="n", ylabel="Length", legend=:none)

This is very chaotic (that's the main feature of the sequence, it is out of control). So we'd better plot the data directly. We can keep the lines as a guide:

plot([collatz(i) for i ∈ 1:100], xlabel="n", ylabel="Length", legend=:none, ls=:dot);
scatter!([collatz(i) for i ∈ 1:100], xlabel="n", ylabel="Length", legend=:none, color=:red)

This is much better than our earlier (text-based) visualization that made us very impressed with 27 (we stopped at 30!) Here we see points with a greater lengths start to become commonplace.

scatter([collatz(i) for i ∈ 1:10000], xlabel="n", ylabel="Length", legend=:none, markersize=.75)

One can also turn to other types of markers:

Plots.supported_markers()
24-element Array{Symbol,1}:
 :none     
 :auto     
 :circle   
 :rect     
 :star5    
 :diamond  
 :hexagon  
 :cross    
 :xcross   
 :utriangle
 :dtriangle
 :rtriangle
 :ltriangle
 :pentagon 
 :heptagon 
 :octagon  
 :star4    
 :star6    
 :star7    
 :star8    
 :vline    
 :hline    
 :+        
 :x 

histogram([collatz(i) for i ∈ 1:250000], bins=0:5:350, norm=true,
    label="Collatz", xlabel="Length", ylabel="Probability")

Let us now look at something which interests us a lot as Physicists: Mathematical functions. Plotting them can be done as follows:

plot(sin, -π, π, lw=5, lc=:blue, legend=:top, label="Sin", xlabel="x")
plot!(x->x, -π, π, label="1st-order")
plot!(x->x-x^3/6, -π, π, label="2nd-order")
plot!(x->x-x^3/6+x^5/120, -π, π, label="3rd-order")

plot(sin, -π, π, lw=5, lc=:blue, legend=:top, dpi=150,
    title="Taylor expansion", label=L"\sin(x)", xlabel=L"x", 
    ylabel=L"\sum_{k=0}^\infty f^{(k)}(0){x^k\over k!}")
plot!(x->x, -π, π, label=L"x")
plot!(x->x-x^3/6, -π, π, label=L"x-x^3/6")
plot!(x->x-x^3/6+x^5/120, -π, π, label=L"x-x^3/6+x^5/120")

plot(sin, x->x, -π/2, π/2, lw=5, xlabel=L"x", ylabel=L"f(x)",dpi=150, label=L"\arcsin")
plot!(sin, x->x, -π, π, legend=:topleft, ls=:dash, color=:blue, label=L"\sin^{-1}", xlim=(-1.1, 1.1), aspectratio=1/2)

Other types of plots:

R(θ) = 3cos(θ)*sin(θ)
plot(R, 0:.01:2π, proj=:polar, lims=(0,1.5), dpi=150, label=L"P_2^1(\theta)")

x=-1:0.01:1
y=-1:0.01:1
h(x,y)=x^2-y^2;

p1=surface(x,y,h,
   xlabel=L"x", ylabel=L"y", zlabel=L"x^2-y^2");
p2=wireframe(x,y,h,
    xlabel=L"x", ylabel=L"y", zlabel=L"x^2-y^2");
plot(p1,p2)

The layout can be changed with layout=(2,1).

A 2D density plot:

f(x,y) = sin(3x) * cos(x+y)

x = 0:0.01:5
y = 0:0.01:5
contour(x, y, f, xlabel="x", ylabel="y", fill=true)

A Rotational field:

xs = -2:0.2:2
ys = -2:0.2:2

df(x, y) = [-y, x] ./ 7

xxs = [x for x in xs for y in ys]
yys = [y for x in xs for y in ys]

quiver(xxs, yys, quiver=df, 
    aspectratio=1,xlims=(-1.75,1.75),ylims=(-1.75,1.75),xlabel=L"x", ylabel=L"y")

Links

• Plots
• https://www.math.purdue.edu/~allen450/Plotting-Tutorial.html
• https://github.com/nassarhuda/JuliaTutorials/blob/master/plotting.ipynb
• Gallery of plots with Plots.jl.