Posts Tagged ‘plot’

gnuplot for dummies…

mercoledì, Settembre 22nd, 2021

At first glance, it might seem that Gnuplot is more difficult to use for plotting than MS Excel. But it only seems that way, the entry threshold is a little higher (you can’t use a mouse to click it, you need to read the documentation here), but in practice, it comes out much easier and more convenient. I wrote a script once and have been using it all my life. It’s really much more difficult for me to plot a graph in Excel, where everything is not logical than in Gnuplot. And the main advantage of Gnuplot is that you can embed it into your programs and visualize data on the fly. Also Gnuplot without any problems builds a graph from a 30-gigabyte file of statistical data, while Exel simply crashed and could not open it.

The pluses of Gnuplot include the fact that it can be easily integrated into the code in standard programming languages. There are ready-made libraries for many languages, I personally came across Fortran and Python. Thus, you can generate graphs directly from your program.

Gnuplot — real application

There are two ways to work with Gnuplot: command mode and script execution mode. I recommend using the second mode immediately, as it is the most convenient and fastest. Moreover, the functionality is absolutely the same.

The graph is plotted by the plot command. As command parameters, you can specify a function or the name of the data file. As well as which column of data to use and how to connect the points, how to denote them, etc. Let’s illustrate.

gnuplot> plot sin(x)
Gnuplot interactive plot of sin(x) function

After executing the command, a window will open where there will be a sine graph with default signatures. Let’s improve this graph and figure out some additional options at the same time.
Let’s sign the axes.

set xlabel "X"

Specifies the label for the abscissa axis.

set ylabel "Y"

Specifies the label for the ordinate axis.

Let’s add a grid to show where the graph is built.

set grid

We will set the limits of the values to which the graph will be limited via the command

set yrange [-1.1:1.1]

Thus, the graph will be drawn from the minimum value of -1.1 to a maximum value of 1.1.

Likewise, I set the range for the abscissa so that only one period of the sinusoid is visible.

set xrange[-pi:pi]

We need to add a title to our schedule:

set title "Gnuplot's examples" font "Helvetica Bold, 20"

Note that you can set fonts and their sizes. See the Gnuplot documentation for what fonts can be used.

And finally, in addition to the sine on the graph, let’s also draw the cosine, and also set the line type and its color. And also add legends, what are we drawing.

plot sin(x) title "sin(x)" lc rgb "red", cos(x) title "cos(x)" lc rgb "green"

Here we draw two graphs on one canvas, in red and green. In general, there are a lot of options for lines (dotted line, stroke, solid), line widths, colors. Like the types of points. My goal is only to demonstrate the range of possibilities. Let’s bring all the commands into one pile and execute them sequentially.

the result:

If you honestly repeated all this after me, you might have noticed that manually entering it every time, even copying, is somehow not comme il faut. But is this a ready-made script? So let’s do it!

a gnuplot script

We make it executable and run it.

chmod +x gnuplottest.gpi

As a result, we get the same window with charts. If you don’t add -persist to the title, the window will automatically close after the script is executed.

In order to output data to a file and not to the screen, you need to reassign the terminal.

set terminal png size 800, 600
set output "result.png"
rep

As you might guess, we indicate the file type, its resolution, then we indicate the file name. Add these two lines to the beginning of our script, and we get this picture in the current folder.

Our final Gnuplot graph

In order to save to postscript, you need to use the following commands:

set terminal postscript eps enhanced color solid
set output "result.ps"
rep

In Python?

Here a simple code for plot a cosine function in Python using matplotlib:

import numpy as np
import matplotlib.pyplot as plt
t = np.arange(0.0,2.0,0.01)
s = 1+ np.cos(2*np.pi*t)
plt.plot(t,s,'*')
plt.xlabel('Time(t)')
plt.ylabel('Voltage(mV)')
plt.title('Cosine wave plot(cos(x))')
plt.show()

Basic Plotting with matplotlib

sabato, Settembre 18th, 2021

The Python ecosystem gives us the option of visualizing relationships between numbers via matplotlib, a plotting library (i.e., not part of core Python) which can produce quality figures. Inside the matplotlib package is the matplotlib.pyplot module, which is used to produce figures in a MATLAB-like environment.

Here a simple example code for the matplotlib.pyplot module.

import matplotlib.pyplot as plt

def plotex(cxs,cys,dxs,dys):
    plt.xlabel('$x$', fontsize=10)
    plt.ylabel('$f(x)$', fontsize=10)
    plt.plot(cxs, cys, 'r-', label='quadratic: $y=x^2$')
    plt.plot(dxs, dys, 'b--^', label='other function: $y=x^{1.8}-1/2$')
    plt.legend()
    plt.show()

x1 = [0.1*i for i in range(60)]
y1 = [x**2 for x in x1]

x2 = [i for i in range(7)]
y2 = [x**1.8 - 0.5 for x in x2]

plotex(x1, y1, x2, y2)
A matplotlib.pyplot module example

We start by importing matplotlib.pyplot in the (standard) way which allows us to use it below without repeated typing of unnecessary characters.

We then define a function, plotex(), that takes care of the plotting, whereas the main program simply introduces four list comprehensions and then calls our function. If you’re still a beginner, you may be wondering why we defined a Python function in this code. An important design principle in computer science goes by the name of separation of concerns (or sometimes information hiding or encapsulation): each aspect of the program should be handled separately. In our case, this means that each component of our task should be handled in a separate function.

Let’s discuss this function in more detail. Its parameters are (meant to be) four lists, namely two pairs of \(x_i\) and \(y_i\) values. The function body starts by using xlabel() and ylabel() to provide labels for the x and y axes. It then creates individual curves/sets of points by using matplotlib’s function plot(), passing in the x-axis values as the first argument and the y-axis values as the second argument. The third positional argument to plot() is the format string: this corresponds to the color and point/line type. In the first case, we used r for red and - for a solid line.

Some of the most important line styles/markers in matplotlib

The fourth argument to plot() is a keyword argument containing the label corresponding to the curve. In the second call to plot() we pass in a different format string and label (and, obviously, different lists); observe that we used two style options in the format string: -- to denote a dashed line and ˆ to denote the points with a triangle marker. The function concludes by calling legend(), which is responsible for making the legend appear, and show(), which makes the plot actually appear on our screen.

We could fine-tune almost all aspects of our plots, including basic things like line width, font size, and so on. For example, we could get LaTeX-like equations by putting dollar signs inside our string, e.g., ‘$x i$’ appears as \(x_i\).

Electric Field of a Distribution of Point Charges

Very briefly, let us recall Coulomb’s law: the force on a test charge Q located at point P (at the position r), coming from a single point charge q0 located at r0 is given by:

\(\mathbf{F}_{0}=k \frac{q_{0} Q}{\left(\mathbf{r}-\mathbf{r}_{0}\right)^{2}} \frac{\mathbf{r}-\mathbf{r}_{0}}{\left|\mathbf{r}-\mathbf{r}_{0}\right|}\)

where Coulomb’s constant is \(k=1 /\left(4 \pi \epsilon_{0}\right)\) in SI units (and \(\epsilon_0\) is the permittivity of free space).

The force is proportional to the product of the two charges, inversely proportional to the square of the distance between the two charges, and points along the line from charge q0 to charge Q. The electric field is then the ratio of the force F0 with the test charge Q in the limit where the magnitude of the test charge goes to zero. In practice, this given us:

\(\mathbf{E}_{0}(\mathbf{r})=k q_{0} \frac{\mathbf{r}-\mathbf{r}_{0}}{\left|\mathbf{r}-\mathbf{r}_{0}\right|^{3}}\)

where we cancelled out the Q and also took the opportunity to combine the two denominators. This is the electric field at the location r due to the point charge q0 at r0.

If we were faced with more than one point charge, we could apply the principle of superposition: the total force on Q is made up of the vector sum of the individual forces acting on Q. As a result, if we were dealing with the n point charges q0, q1,…, qn −1 located at r0, r1,…, rn−1 (respectively) then the electric field at the location r is:

\(\mathbf{E}(\mathbf{r})=\sum_{i=0}^{n-1} \mathbf{E}_{i}(\mathbf{r})=\sum_{i=0}^{n-1} k q_{i} \frac{\mathbf{r}-\mathbf{r}_{i}}{\left|\mathbf{r}-\mathbf{r}_{i}\right|^{3}}\)

Note that you can consider this total electric field at any point in space, r. Note, also, that the electric field is a vector quantity: at any point in space this E has a magnitude and a direction. One way of visualizing vector fields consists of drawing field lines, namely imaginary curves that help us keep track of the direction of the field. More specifically, the tangent of a field line at a given point gives us the direction of the electric field at that point. Field lines do not cross; they start at positive charges (“sources”) and end at negative charges (“sinks”).

We will plot the electric field lines in Python; while more sophisticated ways of visualizing a vector field exist (e.g., line integral convolution), what we describe below should be enough to give you a qualitative feel for things.

We are faced with two tasks: first, we need to produce the electric field (vector) at several points near the charges and, second, we need to plot the field lines in such a way that we can physically interpret what is happening.

Code below is a Python implementation, where Coulomb’s constant is divided out for simplicity.

import numpy as np
import matplotlib.pyplot as plt
from math import sqrt
from copy import deepcopy

def makefield(xs, ys):
    qtopos = {1: (-1,0), -1: (1,0)}
    n = len(xs)
    Exs = [[0. for k in range(n)] for j in range(n)]
    Eys = deepcopy(Exs)
    for j,x in enumerate(xs):
        for k,y in enumerate(ys):
            for q,pos in qtopos.items():
                posx, posy = pos
                R = sqrt((x - posx)**2 + (y - posy)**2)
                Exs[k][j] += q*(x - posx)/R**3
                Eys[k][j] += q*(y - posy)/R**3
    return Exs, Eys

def plotfield(boxl,n):
    xs = [-boxl + i*2*boxl/(n-1) for i in range(n)]
    ys = xs[:]
    Exs, Eys = makefield(xs, ys)
    xs=np.array(xs); ys=np.array(ys)
    Exs=np.array(Exs); Eys=np.array(Eys)
    plt.streamplot(xs, ys, Exs, Eys, density=1.5, color='m')
    plt.xlabel('$x$')
    plt.ylabel('$y$')
    plt.show()

plotfield(2.,20)

We start by importing numpy and matplotlib, since the heavy lifting will be done by the function streamplot(), which expects NumPy arrays as input. We also import the square root and the deepcopy() function, which can create a distinct list-of-lists.

The function makefield() takes in two lists, xs and ys, corresponding to the coordinates at which we wish to evaluate the electric field (x and y together make up r). We also need some way of storing the ri at which the point charges are located. We have opted to store these in a dictionary, which maps from charge qi to position ri. For each position r we need to evaluate E(r): in two dimensions, this is made up of Ex (r) and Ey (r), namely the two Cartesian components of the total electric field. Focusing on only one of these for the moment, say Ex (r), we realize that we need to store its value for any possible r, i.e., for any possible x and y values. We decide to use a list-of-lists, produced by a nested list comprehension. We then create another list-of-lists, for Ey(r). We need to map out (i.e., store) the value of the x and y components of the total electric field, at all the desired values of the vector r, namely, on a two-dimensional grid made up of xs and ys. This entails computing the electric field (contribution from a given point charge qi ) at all possible y’s for a given x, and then iterating over all possible x’s. We also need to iterate over our point charges qi and their locations ri; we do this by saying for q, pos in qtopos.items(): at which point we unpack pos into posx and posy. We thus end up with three nested loops: one over possible x values, one over possible y values, and one over i. All three of these are written idiomatically, employing items() and enumerate().

Our second function, plotfield(), is where we build our two-dimensional grid for the xs and ys. We take in as parameters the length L and the number of points n we wish to use in each dimension and create our xs using a list comprehension; all we’re doing is picking x’s from −L to L. We then create a copy of xs and name it ys. After this,
we call our very own makefield() to produce the two lists-of-lists containing Ex (r) and Ey (r) for many different choices of r.

The result of running this code is shown in the following figure.

Visualizing the electric fields resulting from two point charges.

Basic data manipulation in Python…

martedì, Settembre 14th, 2021

In this post, we will deal with data from ECDC and we will explain basic data manipulation in Python with the Pandas package.

In our day, data is everywhere in enormous size and depth. Data science is an emerging field that penetrates every aspect of our life and, lately, it has proved to be an extraordinary weapon for predicting infections from Covid-19 and organizing strategies to limit the damage.

To import Pandas and Matplotlib packages we code:

import pandas as pd
import matplotlib.pyplot as plt

We download the excel file locally from ECDC site and open it using the read_excel function of Pandas library. We have named the file as data.xls in our case.

df=pd.read_excel("data.xlsx", engine="openpyxl")

We can first explore the data and the columns of the dataframe df:

We observe the columns of the dataframe — in our case, we will use the columns: dateRep, cases and deaths. Additionally, the name of the country is stored in column countriesAndTerritories.

We next select ‘Italy’ as the country under study. A new column is created named DateTime of type datetime where we store the day. In the following, we create a new dataframe with the name df_italia which is the same as the dataframe df_italia_sorted but it is sorted according to the column DateTime

df_italia=df[df.countriesAndTerritories=='Italy']
df_italia['DateTime']=pd.to_datetime(df_italia['dateRep'],format="%d/%m/%Y")

#We sort according to DateTime
df_italia_sorted=df_italia.sort_values(by='DateTime')

df_italia_selected=df_italia_sorted[df_italia_sorted.month>4]

We are interested in data after the month of April (i.e., May, June, July, August, … etc) so we choose to filter using the column month and create a new dataframe df_italia_selected.

Since the data in columns cases and deaths may have great variation, it is practical in order to understand the trend to use a moving average. We choose a moving average of seven days and we create two new columns (Moving Average Cases and Moving Average Deaths) where we store the average values of cases and deaths.

#Calculate moving average

df_italia_selected['Moving Average Cases']=df_italia_selected.cases.rolling(7,min_periods=1).mean()
df_italia_selected['Moving Average Deaths']=df_italia_selected.deaths.rolling(7,min_periods=1).mean()

We now plot the cases and deaths as functions of time. We choose the red color for cases and blue for deaths. It is useful to plot cases and deaths in the same figure with common x-axis in order to understand possible connection and relation. So, we use the subplots function and first create figure fig and axis ax1 (this will be the axis for the cases and it will be the left axis). We then create ax2 using twinx function. The values for deaths will be our right axis. A dashed line is used for the average values.

#Figure

fig, ax1=plt.subplots()
color1='tab:red'

ax1.plot(df_italia_selected['DateTime'], df_italia_selected['cases'], color=color1)

ax1.plot(df_italia_selected['DateTime'], df_italia_selected['Moving Average Cases'], color=color1,linestyle='dashed')

ax1.set_xlabel('Data')

ax1.set_ylabel('Cases',color=color1)
ax1.tick_params(axis='y',labelcolor=color1)

locs, labels=plt.xticks()
plt.setp(labels,rotation=90)

ax2=ax1.twinx() #instantiate a second axes that shares the same x-axis
color2='tab:blue'

ax2.plot(df_italia_selected['DateTime'], df_italia_selected['deaths'], color=color2)
ax2.plot(df_italia_selected['DateTime'], df_italia_selected['Moving Average Deaths'], color=color2,linestyle='dashed')

ax2.set_ylabel('Deaths',color=color2)
ax2.tick_params(axis='y',labelcolor=color2)

fig.tight_layout() #otherwise the right y-label is slightly clipped

The figure below is the program output.

Cases and deaths as a function of data for Italy

Python animations with Matplotlib

mercoledì, Settembre 1st, 2021

In Python, plotting graphs is straightforward — you can use powerful libraries like Matplotlib. It happens, however, that you need to visualize the trend over time of some variables – that is, you need to animate the graphs.

Luckily, it’s just as easy to create animations as it is to create plots with Matplotlib.

Matplotlib

Matplotlib – as you can read on the official site – is a comprehensive library for creating static, animated, and interactive visualizations in Python. You can plot interactive graphs, histograms, bar charts, and so on.

How to Install Matplotlib

Installing Matplotlib is simple. Just open up your terminal and run:

pip install matplotlib

Numpy

Also, if you don’t have numpy, please install it so you can follow the examples in this tutorial:

pip install numpy

How to Plot with Matplotlib

Even though this tutorial is about animations in Matplotlib, first let’s create a simple static graph of a sine wave:

import numpy as np
import matplotlib.pyplot as plt

x = np.arange(0, 10, 0.1)
y = np.sin(x)

fig = plt.figure()
ax = plt.axes(xlim=(0, 10), ylim=(-1.1, 1.1))
diagram = plt.plot(x, y)
 
plt.show()
A basic sine wave

How to Animate with Matplotlib

To create an animation with Matplotlib you need to use Matplotlib’s animation framework’s FuncAnimation class.

For instance, let’s create an animation of a sine wave:

import numpy as np
from matplotlib import pyplot as plt
from matplotlib.animation import FuncAnimation

fig = plt.figure()
ax = plt.axes(xlim=(0, 5), ylim=(-1.5, 1.5))
line, = ax.plot([], [], lw=2)

def init():
    line.set_data([], [])
    return line,

def animate(i):
    x = np.linspace(0, 5, 500)
    y = np.sin(2 * np.pi * (x + 0.02 * i))
    line.set_data(x, y)
    return line,

sinewawe_animation = FuncAnimation(fig, animate, init_func=init, frames=200, interval=20, blit=True)
sinewawe_animation.save("Animation.gif")

plt.show()

We have:

Let’s then go through the code above in a bit more detail to better understand how animations work with Matplotlib.

Lines 1–3

Here you add the required libraries. In particular, we add the FuncAnimation class that can be used to create an animation for you.

Lines 5–7

fig = plt.figure()
ax = plt.axes(xlim=(0, 5), ylim=(-1.5, 1.5))
line, = ax.plot([], [], lw=2)

Here you first create an empty window for the animation figure. Then you create an empty line object. This line is later modified to form the animation.

Lines 9–11

def init():
    line.set_data([], [])
    return line,

Here you create an init() function that sets the initial state for the animation.

Lines 13–17

You then create an animate() function. This is the function that gives rise to the sine wave. It takes the frame number i as its argument, then it creates a sine wave that is shifted according to the frame number (the bigger it is, the more the wave is shifted). Finally, it returns the updated line object. Now the animation framework updates the graph based on how the line has changed.

def animate(i):
    x = np.linspace(0, 5, 500)
    y = np.sin(2 * np.pi * (x + 0.02 * i))
    line.set_data(x, y)
    return line,

Line 19

sinewawe_animation = FuncAnimation(fig, animate, init_func=init, frames=200, interval=20, blit=True)

This line of code puts it all together and creates the actual animation object. It simply:

  • Renders an animation to the figure fig by repeatedly calling the animate() function starting from the initial state defined by init()
  • The number of frames rendered to “one round of animation” is 200.
  • A delay between two frames is 20 milliseconds (1000ms / 20ms = 50 FPS).
  • (The blit=True makes sure only changed pieces of the plot are drawn, to improve the efficiency)

Line 21

sinewawe_animation.save("Animation.gif")

This piece of code is used to generate an animated gif (the same one I used in this tutorial to show you the animation)

Line 22

You guessed it — this line shows the animation.