Gnuplot tricks

Projecting contours

2010-09-21T14:17:00.000-07:00

Karl asked a question some time ago, in which he wanted to know how one can produce this graph. As I pointed out in my reply, it is rather easy, if we can rotate the data file by 90 degrees. I will only post a skeleton here, you can dress up the graph at your will.

For a start, here is our data file, which we will call 'out.dat'

-0.299  -0.265  -0.215  -0.151  -0.078  0.000   0.078   0.151   0.215   0.265   0.299
-0.513  -0.455  -0.368  -0.259  -0.134  0.000   0.134   0.259   0.368   0.455   0.513
-0.694  -0.616  -0.499  -0.351  -0.181  0.000   0.181   0.351   0.499   0.616   0.694
-0.833  -0.738  -0.596  -0.411  -0.191  0.037   0.243   0.430   0.600   0.739   0.833
-0.919  -0.812  -0.624  -0.271  0.287   0.736   0.767   0.658   0.697   0.819   0.920
-0.949  -0.832  -0.582  0.048   1.186   2.000   1.680   1.007   0.781   0.851   0.949
-0.919  -0.812  -0.624  -0.271  0.287   0.736   0.767   0.658   0.697   0.819   0.920
-0.833  -0.738  -0.596  -0.411  -0.191  0.037   0.243   0.430   0.600   0.739   0.833
-0.694  -0.616  -0.499  -0.351  -0.181  0.000   0.181   0.351   0.499   0.616   0.694
-0.513  -0.455  -0.368  -0.259  -0.134  0.000   0.134   0.259   0.368   0.455   0.513
-0.299  -0.265  -0.215  -0.151  -0.078  0.000   0.078   0.151   0.215   0.265   0.299

and its "rotated" pair, 'out2.dat'

-0.299  -0.513  -0.694  -0.833  -0.919  -0.949  -0.919  -0.833  -0.694  -0.513  -0.299
-0.265  -0.455  -0.616  -0.738  -0.812  -0.832  -0.812  -0.738  -0.616  -0.455  -0.265
-0.215  -0.368  -0.499  -0.596  -0.624  -0.582  -0.624  -0.596  -0.499  -0.368  -0.215
-0.151  -0.259  -0.351  -0.411  -0.271  0.048   -0.271  -0.411  -0.351  -0.259  -0.151
-0.078  -0.134  -0.181  -0.191  0.287   1.186   0.287   -0.191  -0.181  -0.134  -0.078
0.000   0.000   0.000   0.037   0.736   2.000   0.736   0.037   0.000   0.000   0.000
0.078   0.134   0.181   0.243   0.767   1.680   0.767   0.243   0.181   0.134   0.078
0.151   0.259   0.351   0.430   0.658   1.007   0.658   0.430   0.351   0.259   0.151
0.215   0.368   0.499   0.600   0.697   0.781   0.697   0.600   0.499   0.368   0.215
0.265   0.455   0.616   0.739   0.819   0.851   0.819   0.739   0.616   0.455   0.265
0.299   0.513   0.694   0.833   0.920   0.949   0.920   0.833   0.694   0.513   0.299

These were produced in octave by the function

f(x,y) = sin(y/4)*cos(x/4)+exp(-x*x - y*y/3)

Instead of actually rotating the date file, I simply interchanged the variables, and printed out the file for a second time, for I was a bit lazy...

Anyway, this is what we have to do:

reset
unset key
set contour base
set pm3d at ss

set xrange [-2:10]
set yrange [0:12]

splot for [i=1:10:2] 'out.dat' u (-2):0:i w l lt i, \
for [i=1:10:2] 'out2.dat' u 0:(12):i w l lt i,\
'out.dat' matrix w pm3d

This is really simple: We plot the contours by plotting 'out.dat', and 'out2.dat' column by column, and keeping the first and second coordinates constant. In this way, we "project" those columns onto the y-z and x-z planes. In order to make the contours more visible, we have to specify an xrange and yrange which is a bit bigger, than our actual data set. At the end, we plot the data file as surface. If we set the contour beforehand, we will see the contours on the bottom.
And here is the figure that we have just produced. Don't be fooled by the fact that there are only three lines on the x-z plane: since our function was symmetric in y with respect to, some contour lines will overlap. And again, this graph should still be properly annotated.

A small (or big) diversion

2010-08-26T00:34:00.000-07:00

In the past year, I have been trying to argue on these pages that gnuplot has some advantages over many other plotting utilities. Ultimate control over graph properties, the simplicity of plotting, the ease of scripting. These are the strengths of gnuplot that spring to mind first. At least, to my mind. At the same time, I also have to admit that there are weaknesses. And these weaknesses all come down to the same deficiency: the lack of a certain modularity, both at the user level, and in the code. This makes expanding gnuplot extremely difficult, at times, impossible. At the user level, we have to use what we have, and there are not too many options when it comes to even such simply tasks as calculating the average of a data set. If something is not implemented in the code, the user is not "supposed" to use it. Now, in gnuplot 4.4, some of these problems can be overcome with some witty scripting, and mainly, abuse of procedures. If you want to see some nasty hacks, just skim through these pages. And these issues all become even more problematic, when it comes to trying to fix the problem at the developer's level. The code, as it is written now, does not support straightforward expansion, even implementing as simple things as, again, calculating the average of a data set are somewhat tricky. At yet another level, even if the code is fixed, the original gnuplot code is not published under GPL, therefore, changing it does not mean that the battle is won: one can't just take the code, modify it, and put it up on a web page.

These were the problems that I have realised in the past couple of years, and this is why I decided to re-think certain things, and start the development of a plotting utility. I wanted to keep what was good in gnuplot, but I wanted to right the wrongs. I was seriously pondering where and how to set out, when someone pointed out to me that I am not the first person, who faces the same dilemma, and that there is another project already underway. In fact, the other project was quite advanced, when I caught glimpse of it. And I have to say that what I saw was rather impressive. It impressed me, because its development is done along the lines that I mentioned above, it is thought-over well, and it is already quite mature. It can easily compete with gnuplot, for most things are already implemented, and it has the modularity that I missed so much. And it is under GPL, so you can do whatever you want. Well, almost.

With these remarks, I would like to call your attention to pyxplot. If you haven't seen it yet, please, visit the web site, and give it a try! Their main site is here, and you can find a number of very pleasing plots here. Most of your gnuplot scripts will work "out of the box", and those that won't, can be tweaked very easily. At the end of the pdf manual, you can find a discussion on what the differences between gnuplot and pyxplot are, and how you can make your scripts work. At the same time, enjoy the convenience of easy unit manipulations, data analysis, things like Fourier transforms and filtering, numerical integration and differentiation, the option for defining not only simple functions but procedures for your common tasks.

As a concluding remark, I would also like to announce a parallel blog of mine, pyxplot-tricks, where I will discuss how and what can be achieved in pyxplot. I will still try to keep gnuplot-tricks active, and I will certainly answer questions posted here.
Cheers,
Zoltán

Fence plots with a some-liner

2010-07-14T03:58:00.000-07:00

About this time last year, I showed how one can produce fence plots in gnuplot, even if the data is from a file, not from a function. (The function plot is somewhat trivial, you can find it amongst the demos.) I used some heavy data processing then, though everything was handled in gnuplot. With the arrival of version 4.4, all that machinery can be made much simpler. If you continue on reading, you fill find a quite straightforward and flexible method.

For a start, we will need some data. Instead of generating it in gnuplot, I will just post my data file, which reads as follows

1 2 3 1 2 5
2 2 4 2 3 1
3 4 3 6 1 1
4 2 3 1 1 4
5 3 2 5 4 3
6 2 3 6 5 3

We have six columns here, but only five are the data: the first columns is for indexing, or whatever you like. This does not change the idea. In what follows, I will call this file '3fill.dat'

Now, our first script looks like this

reset
unset key
unset colorbox
set ytics offset 0,-1
set ticslevel 0
min = 0
col = 5

DATA = ""
DATA2 = ""
PALETTE = "set palette defined ("

pr(x, y) = sprintf("%f %f\n", x, y)
zero_line(x, y) = DATA.sprintf("\n").DATA2.sprintf("\n%f %f\n", x, y)
zero_pal(x) = sprintf("%d %.3f %.3f %.3f", x, rand(0), rand(0), rand(0))

f(x, y) = ($0 == 0 ? (DATA = zero_line($1, x), DATA2 = pr($1, min), PALETTE = PALETTE.zero_pal(y).", ") : \
        (DATA = DATA.pr($1, x), DATA2 = DATA2.pr($1, min)), x)

plot for [i=2:col+1] '3fill.dat' u 1:(f(column(i), i))

DATA = DATA.sprintf("\n").DATA2

set print '3fill.tab'
print DATA
set print

eval(PALETTE.zero_pal(col+2).")")

splot for [i=0:col-1] '3fill.tab' every :::(2*i)::(2*i+1) u 1:(i):2:(i+2) w pm3d

and the figure produced is here

Once you absorb it, the script is really simple. The first line where something actually happens is where we define DATA, DATA2, and PALETTE. What we will do is to read in the data from the file, and then add the numbers to a string. Once we have read all data, we print the string to a file, and then use that file as our new data file. The reason for this is that we have, in some sense, for to duplicate our data: in each column, we have one set of data, and what this set determines is a curve. In order to produce a fence, however, we need a surface. The simples way of getting this surface is to print the data twice. Of course, we have to modify the data a bit, but this is the only trick here.

We define three functions, pr, which is just a short-hand for formatted printing of two numbers, zero_line, which is again, a printing routine, when the record number is zero, i.e., when we are processing the first data point in each column, and zero_pal, which generates a new palette colour, as we enter a new column. If you are satisfied with some readily available palette, you can skip this function, and any calls to it.

Next, we define a function, f(x,y), which makes use of the three above-mentioned functions, and amounts to the data-duplication process. In order to reduce the complexity of the problem, we will print each column, and its duplicate in the same file. This, however, means that we have got to differentiate between various data sets. We do this by inserting and extra blank line each time we are faced with a new column. This is why we have to distinguish the $0 == 0 case in f(x,y). Also, the palette has to be re-defined only if there is a new data set.

If you watch carefully here, there are two strings, DATA, and DATA2. DATA2 is to hold the duplicate, while DATA is an ever-expanding string with the original, and the duplicate data. In the duplicate, we don't actually hold any duplicate, we simply print the indices (these are in the first column), and a constant number, min, which we defined at the very beginning. The value of this determines how tall (or how deep) the fences will be.

In order to generate the duplicated data set, first we call a dummy plot with f(x,y), add the last duplicate, DATA2 to DATA (this is required, because we concatenate DATA and DATA2 only, when $0 == 0, i.e., the last data set is not added to DATA automatically. Having defined DATA, we print the everything to a file. Note that we could process a number of columns easily, thanks to the for loop in the dummy plot.

At this point, we have everything in the new data file, we have only got to plot it. Recall, that all columns are in one file now, and we have to separate them in the new plot. This is why we use the 'every' keyword when stepping through the data sets.

We can easily add a grid to the figure, if we call another plot at the very end of our script: just add

for [i=0:col-1] '3fill.tab' every :::(2*i)::(2*i+1) u 1:(i):2 w l lt -1

after the last line, and get the following figure

If the order of the plots is exchanged, the "wires" will be covered by the panes of the fence, so it will not give this impression of semi-transparency. This is it for today. In my next post, I will elaborate on for what else we can use the trick above.
Cheers,
Zoltán

Broken axis, once more

2010-06-13T03:55:00.000-07:00

I have discussed this subject at least on two occasions, and in fact, most of the present post was already described in one of my very first posts. However, I thought that it might be worthwhile to dust it off, especially, that this is a scientifically relevant feature, which is still missing in gnuplot. But with a bit of work, we can rectify the problem. In short, there are situations, when we just have to break the axis, simply because there is a large gap between two relevant ranges of a plot. Take this example, for instance

plot [0:40] 20.0*atan(x-20.0) + 32 + sin(x)

which would result in the following graph:

The problem is obvious: the function has some interesting modulation close to 0, and close to 60, but it is rather dull between these two extrema. The solution is to cut out the segment between 10, and 60, say.

What we will use is a very handy function in gnuplot 4.4, which lets the user set the position of the graph exactly. In a multiplot, we would usually set the position as

set multiplot 
set size 1, 0.3
set origin 0, 0.5
plot 'foo' using 1:2 

set size 1, 0.2
set origin 0, 0
plot 'bar' using 3:4

which produces two graphs of size 1, 0.3, and 1, 0.2, respectively, and places them in such a way that their bottom left corner is at (0, 0.5), and (0, 0). But when we say "the bottom left corner", we actually mean the whole figure, tic marks, axis labels, everything. This means at least two things. One is that if we want to break the axis using a multiplot, the ranges will not necessarily be proportional on the figure, simply because the size referred to the size of the whole graph, and not only to the plotting area. Second, if the size of the tic labels is different in the two graphs, they will no longer be aligned properly. This would happen, e.g., if we were to plot over the ranges [0:9], and [1000:10009]: the first range requires labels of width 1, while the second labels of width 5, therefore, the second graph would be narrower, and its left vertical axis shifted to the right, at least, as far as the plotting area is concerned.

In gnuplot 4.4, however, one can set the positions of the plots, and not the whole graph. This is achieved by issuing a command similar to this

set lmargin at screen 0.1

which aligns the left vertical axis of the graph with the point that is at screen position 0.1 along the horizontal direction. There are three more margins, rmargin, tmargin, and bmargin, setting the right hand side, the top, and the bottom of the graph. Specifying the plot's corners explicitly removes the above-mentioned problem with the alignments.

Having said this, our script could read as follows

reset
unset key
bm = 0.15
lm = 0.12
rm = 0.95
gap = 0.03
size = 0.75
y1 = 0.0; y2 = 11.5; y3 = 58.5; y4 = 64.0

set multiplot
set xlabel 'Time [ns]'
set border 1+2+8
set xtics nomirror
set ytics nomirror
set lmargin at screen lm
set rmargin at screen rm
set bmargin at screen bm
set tmargin at screen bm + size * (abs(y2-y1) / (abs(y2-y1) + abs(y4-y3) ) )

set yrange [y1:y2]
plot [0:40] 20.0*atan(x-20.0) + 32 + sin(x)

unset xtics
unset xlabel
set border 2+4+8
set bmargin at screen bm + size * (abs(y2-y1) / (abs(y2-y1) + abs(y4-y3) ) ) + gap
set tmargin at screen bm + size + gap
set yrange [y3:y4]

set label 'Power [mW]' at screen 0.03, bm + 0.5 * (size + gap) offset 0,-strlen("Power [mW]")/4.0 rotate by 90

set arrow from screen lm - gap / 4.0, bm + size * (abs(y2-y1) / (abs(y2-y1)+abs(y4-y3) ) ) - gap / 4.0 to screen \
lm + gap / 4.0, bm + size * (abs(y2-y1) / (abs(y2-y1) + abs(y4-y3) ) ) + gap / 4.0 nohead

set arrow from screen lm - gap / 4.0, bm + size * (abs(y2-y1) / (abs(y2-y1)+abs(y4-y3) ) ) - gap / 4.0  + gap to screen \
lm + gap / 4.0, bm + size * (abs(y2-y1) / (abs(y2-y1) + abs(y4-y3) ) ) + gap / 4.0 + gap nohead

set arrow from screen rm - gap / 4.0, bm + size * (abs(y2-y1) / (abs(y2-y1)+abs(y4-y3) ) ) - gap / 4.0 to screen \
rm + gap / 4.0, bm + size * (abs(y2-y1) / (abs(y2-y1) + abs(y4-y3) ) ) + gap / 4.0 nohead

set arrow from screen rm - gap / 4.0, bm + size * (abs(y2-y1) / (abs(y2-y1)+abs(y4-y3) ) ) - gap / 4.0  + gap to screen \
rm + gap / 4.0, bm + size * (abs(y2-y1) / (abs(y2-y1) + abs(y4-y3) ) ) + gap / 4.0 + gap nohead

plot [0:40] 20.0*atan(x-20.0) + 32 + sin(x)

unset multiplot

The first couple of lines specify how big a figure we want to have: bm, lm, and rm are the bottom, left, and right margins, respectively. We also define the size of the gap, which we will have between the two plots. y1 through y4 are the definitions of our plot ranges. In other words, the interval between y2, and y3 will be cut out of our figure.

In the multiplot environment, we set the axes (for the bottom figure on the bottom, left, and right, while for the top figure on the top, left, and right), the axis labels (the vertical label we have to set by hand, for otherwise it would be centred on the vertical axis of the bottom or top figure, but not on the whole), and set the positions of the figures. Note that the definition used for tmargin and bmargin makes sure that the two plotted intervals are proportional. Before plotting the second curve, we also set four small headless arrows, which are meant to represent the break in the axes. It can be left out, if not desired, or they can be replaced by two dashed vertical lines.

This method can also be used, if one wants to plot a single curve or data set, but with logarithmic axis on one interval, and linear on the other.

Ministry of Sillier Walks

2010-06-12T15:08:00.000-07:00

In my last post, I have shown how we can define or read an array of number from a file. Having constructed the array, the question naturally arises: can we use it for something else, to do something that we could not do otherwise. The answer is yes, although in the present post, there will not be anything that I haven't discussed here or there. We have only got to put the pieces together, and with a little bit of work, we can produce three-dimensional column stacked histograms at ease. I will use the same datafile, but I repeat it here:

"" France Germany Japan Nauru
Defense 9163 4857 2648 9437
Agriculture 3547 5378 1831 1948
Education 7722 7445 731 9822
Industry 4837 147 3449 6111
"Silly walks" 3441 7297 308 7386

Then, let us see, how we are going to make that histogram! For the sake of example, we will draw four cylinders, which will be striped according to the values that the columns in the data file take. We could do two things here. Provided that we have already read the values into an array, we could draw 4 times 5 cylinders of various size and colour. Alternatively, we could draw 4 cylinders, and colour them with stripes that we take from the palette. In this case, we have to find some way of defining the proper palette. Both methods have advantages and disadvantages. The advantage of the first one is that it is faster, because the cylinders are monocolour, which means that we need only two isosamples. On the other hand, we have to have "nested" for loops (which we would just simply write out). The advantage of the second method is that we can get away with one for loop, but we need many isosamples, because the cylinders contain many colour, though the transition between the colours is required to be abrupt. But since we do not know where the boundary between two colours is, we would need many isosamples. Consequently, it will be slow.

When defining our array, we will employ the trick from the last post, and we will do something very similar, when we determine our palette. In addition, we will also have to calculate the sum in the columns, for the cylinders' height should be proportional to that. If, on the other hand, we opt for re-scaling the cylinders to the same height, we, again, need the sum. With these preliminary remarks, the first version of our script could look like this

reset

unset key
unset colorbox
unset xtics; unset xlabel
unset ytics; unset ylabel
unset ztics
file = 'marimekko.dat'   
cylinder = 'cylinder.dat'

set ticslevel 0   
set border 1+2+4+8
set parametric; set urange [0:2*pi]; set vrange [0:1]; set iso 2, 200
set table cylinder
splot cos(u), sin(u), v, cos(u)*v, sin(u)*v, 1
unset table
unset parametric

col = 4 
row = 5 
sm = 0.0
g(x,a) = (abs(x-a) < 0.1 ? 1 : 0)
h(x,a) = (x <= a ? 0.0 : 1.0)

ARRAY = "b(x, y) = 0"
SUM = "s(x,a) = 0"
PALETTE = "set palette defined (-1 0.7 0 0"

array(x, c) = (sm = h($0,0.5)*sm + x, ARRAY = ARRAY.sprintf(" + g(x,%d)*h(y,%.3f)", c, sm), \
                SUM = SUM.sprintf(" + %f*g(x,%d)", x, c), x )
pal(x, c) = (PALETTE = PALETTE.sprintf(", %.1f %.3f %.3f %.3f", c-0.1, rand(0)*0.7, rand(0)*0.7, rand(0)*0.7), x
ff(x, c) = (array(x, c), x)

plot for [i=2:col+1] file every ::1 using 0:( ff(column(i), i) )
plot file every ::1 using 0:(pal(1, column(0)))

eval(ARRAY)
eval(PALETTE.")")
eval(SUM)
         
set xrange [4:5+3*col]
set yrange [-10:3*col-9]

splot for [i=2:col+1] cylinder using ($1+3*i):2:($3*s(i,i)):(b(i,$3*s(i,i))) with pm3d

The first couple of lines set up the graph, and they are trivial, as is the plot to the 'cylinder.dat'. The definition of g(x,a) should be familiar from the last post, and h(x,a) is nothing but the Heaviside function. We have to deal with a matrix, therefore, the string ARRAY begins with "b(x,y) = 0", which will, in due course, become the definition of a two-variable function. SUM, that will bring s(x,a) to life, also appears to define a two-variable function, but this is only apparent: s(x,a) is exactly as much of a two-variable function as is g(x,a). A mere convenience, nothing more. We also begin the definition of a palette, but we stop short of its completion: that will be done during the first plot.

In the array(x,c) function, we read in the numbers from the file, and at the same time, we also calculate the sums in the columns. This function is really similar to that from the last post. We also define a function for filling up the palette. For now, this function writes random colours into the palette at every integer.

Having defined the functions, we "plot" our data file, so as to read in the numbers. We plot the second column of the same file again, in order to prepare our palette. Note that this latter plot is really nothing but a strange way of creating a for loop: we do something as many times as there are elements in a column. In both plots, by applying the every keyword, we skip the first line, which is the column header.

At this point we are almost done: the only remaining thing is the evaluation of our new definitions, and the actual plot. The setting of the xrange and yrange is necessary only in order to ensure that the cylinders are cylinders, not elliptical.

Our first script results in the graph below:

This is OK for a start, but could we improve it a bit? With some work, we could. For one thing, we could add the sum at the top of each cylinder. This can easily be done by augmenting the last plot with the line

for [i=2:col+1] file every ::0::0 using (3*i):(0):(s(i,i)+5e3):(sprintf("%d", s(i,i))) w labels

We can also give a title to each cylinder by reading out the values in the first row. This is done by adding

for [i=2:col+1] file every ::0::0 using (3*i):(-2):(0):(stringcolumn(i)) w labels centre

Finally, we can easily add a legend to the figure. All we have to do is to read out the first column of our data file, and draw 5 cylinders with the appropriate colour. We can achieve this by invoking

file using (0):(0):(2e4-$0*5e3):1 w labels right, \
for [i=1:row] cylinder using ($1+5):($2-5.0):($3*2e3+i*5e3):(i-1) with pm3d

in the last plot. With these modifications, the complete script would look like this

reset

unset key
unset colorbox
unset xtics; unset xlabel
unset ytics; unset ylabel
unset ztics
file = 'marimekko.dat'   
cylinder = 'cylinder.dat'

set ticslevel 0   
set border 1+2+4+8
set parametric; set urange [0:2*pi]; set vrange [0:1]; set iso 2, 200
set table cylinder
splot cos(u), sin(u), v, cos(u)*v, sin(u)*v, 1
unset table
unset parametric

col = 4 
row = 5 
sm = 0.0
g(x,a) = (abs(x-a) < 0.1 ? 1 : 0)
h(x,a) = (x <= a ? 0.0 : 1.0)

ARRAY = "b(x, y) = 0"
SUM = "s(x,a) = 0"
PALETTE = "set palette defined (-1 0.7 0 0"

array(x, c) = (sm = h($0,0.5)*sm + x, ARRAY = ARRAY.sprintf(" + g(x,%d)*h(y,%.3f)", c, sm), \
                SUM = SUM.sprintf(" + %f*g(x,%d)", x, c), x )
pal(x, c) = (PALETTE = PALETTE.sprintf(", %.1f %.3f %.3f %.3f", c-0.1, rand(0)*0.7, rand(0)*0.7, rand(0)*0.7), x
ff(x, c) = (array(x, c), x)

plot for [i=2:col+1] file every ::1 using 0:( ff(column(i), i) )
plot file every ::1 using 0:(pal(1, column(0)))

eval(ARRAY)
eval(PALETTE.")")
eval(SUM)
         
set xrange [4:5+3*col]
set yrange [-10:3*col-9]

splot for [i=2:col+1] cylinder using ($1+3*i):2:($3*s(i,i)):(b(i,$3*s(i,i))) with pm3d, \
for [i=2:col+1] file every ::0::0 using (3*i):(0):(s(i,i)+5e3):(sprintf("%d", s(i,i))) w labels, \              
file using (0):(0):(2e4-$0*5e3):1 w labels right, \
for [i=1:row] cylinder using ($1+5):($2-5.0):($3*2e3+i*5e3):(i-1) with pm3d, \  
for [i=2:col+1] file every ::0::0 using (3*i):(-2):(0):(stringcolumn(i)) w labels centre

and would the following figure:

Well, this is sort of OK, but what if we still do not like it? There are two things that we could easily implement, and would change the character of our graph completely. One is that we can give the cylinders a true 3D lookout, by adding phongs to them. The other one is that we could remove the legends, and add it to one of the cylinders.

So, let us see what we could do in the way of phonging. This is really simple: if we think about it, the phong is nothing but a white spot on our graph, where the saturation of the colours increases towards the centre of the spot. Therefore, all we have to do is to insert white into the palette, but to do it in a way that white saturates all little cylinders. We could, then, modify our palette function as

pal(x, c) = (PALETTE = PALETTE.sprintf(", %.1f %.3f %.3f %.3f, %.1f 1 1 1", \
      c-0.01, rand(0)*0.7, rand(0)*0.7, rand(0)*0.7, c+0.99), x)

and add a function that changes the saturation as

colour(x,y) = 0.25*exp(-(x-0.7)**2/0.2-(y+0.7)**2/0.2)

If we look at the palette function, the random numbers will be between 0-0.7, and the function colour(x,y) will add 0.25 to those numbers at the centre of the spot, which, in this particular case, will be at 45 degrees with respect to the x axis. When plotting, we have to add this function to our cylinders.

As for the labels, we might want to place them at the centre of each coloured cylinder in the last cylinder, representing Nauru. That is, we use the first column, and add the labels from that at half of the height of the cylinders whose size is read from the fifth column. We could use this function

lab(x) = (sm = sm + x, sm-0.5*x)

The value of sm is updated when a new value is read from the fifth column, and the return value of the function is just the cumulative sum minus half of the last value. Note that since we use sm, which was also utilised in array(x,c), we will have to re-set its value to zero before we use it.

With these modifications, the complete script reads as

reset

unset key
unset colorbox
unset xtics; unset xlabel
unset ytics; unset ylabel
unset ztics
file = 'marimekko.dat'
cylinder = 'cylinder.dat'

set ticslevel 0
set border 1+2+4+8
set parametric; set urange [0:2*pi]; set vrange [0:1]; set iso 2, 200
set table cylinder
splot cos(u), sin(u), v, cos(u)*v, sin(u)*v, 1
unset table
unset parametric

col = 4
row = 5
sm = 0.0
g(x,a) = (abs(x-a) < 0.1 ? 1 : 0)
h(x,a) = (x <= a ? 0.0 : 1.0)

ARRAY = "b(x, y) = 0"
SUM = "s(x,a) = 0"
PALETTE = "set palette defined (-1 0.7 0 0, -0.5 1 1 1"

array(x, c) = (sm = h($0,0.5)*sm + x, ARRAY = ARRAY.sprintf(" + g(x,%d)*h(y,%.3f)", c, sm), \
  SUM = SUM.sprintf(" + %f*g(x,%d)", x, c), x )
pal(x, c) = (PALETTE = PALETTE.sprintf(", %.1f %.3f %.3f %.3f, %.1f 1 1 1", \
      c-0.01, rand(0)*0.7, rand(0)*0.7, rand(0)*0.7, c+0.99), x)
ff(x, c) = (array(x, c), x)

colour(x,y) = 0.25*exp(-(x-0.7)**2/0.2-(y+0.7)**2/0.2)
lab(x) = (sm = sm + x, sm-0.5*x)
 
plot for [i=2:col+1] file every ::1 using 0:( ff(column(i), i) )
plot file every ::1 using 0:(pal(1, column(0))) 

eval(ARRAY)
eval(PALETTE.")")
eval(SUM)

set xrange [4:5+3*col]
set yrange [-10:3*col-9]

splot for [i=2:col+1] cylinder using ($1+3*i):2:($3*s(i,i)):(b(i,$3*s(i,i))+colour($1,$2)) with pm3d, \
for [i=2:col+1] file every ::0::0 using (3*i):(0):(s(i,i)+5e3):(sprintf("%d", s(i,i))) w labels, \
sm = 0, file using (3*col+4):(1):(lab($5)):1 w labels left, \
for [i=2:col+1] file every ::0::0 using (3*i):(-2):(0):(stringcolumn(i)) w labels right rotate by 30

I think I cannot add more to this figure, we have explored and exhausted all possibilities here. Next time I will come back to an older topic from this blog, and show how that can be done in an elegant way, applying the new functionalities of gnuplot 4.4.
Cheers,
Zoltán

Ministry of Silly Walks

2010-05-09T11:15:00.000-07:00

In a comment last week, someone asked whether it was possible to draw a Marimekko plot, i.e., a histogram in which both directions contain relevant information. In other words, the question is whether we could draw a square, and populate it with rectangles in such a way that the area of the rectangles is read from a file. I thought that it should be possible, but on the way, I also found out a couple of interesting things. If you keep reading, you will see it for yourself, how we can define a vector, whose value is taken from a data file, and how we can manipulate the elements of that vector. In some sense, this is similar to the trick that we made use of, when we generated a parametric plot from a file. But we can do things in a slightly better fashion.

First, we will need a data file, and for the sake of conformity with the question of the commenter, I will just use this

"" France Germany Japan Nauru
Defense 9163 4857 2648 9437
Agriculture 3547 5378 1831 1948
Education 7722 7445 731 9822
Industry 4837 147 3449 6111
"Silly walks" 3441 7297 308 7386

(We can already deduce that with the sole exception of Japan, countries spend a large chunk of their GDP on silly walk.)

Now, our first attempt could be this:

reset

file = 'marimekko.dat'

set style data histograms
set style histogram columnstacked
set style fill solid border -1
set boxwidth 1.0

set xrange [-1:5]
set yrange [0:5e4]

plot newhistogram at 0, file u 2 title col, \
newhistogram at 1, file u 3 title col, \
newhistogram at 2, file u 4 title col, \
newhistogram at 3, file  u 5 title col, \

and it should be quite obvious that this is not what we want:

It just falls short of our expectations in every respect: There is a gap between the columns, the colours are not consistent, and the width of the columns is equal. The only reasonable thing that happened here is that we can actually set the position of the columns. This will become important later on.

Let us try to improve on the figure, step by step. First, we will place the histograms in a multiplot, for that will make life a lot easier: this is our only way of manipulating the column width during the plot. In this spirit, our second script will be this:

reset

file = 'marimekko.dat'

set style data histograms
set style histogram columnstacked
set style fill solid border -1
set boxwidth 1.0

set xrange [-1:5]
set yrange [0:5e4]
set multiplot

plot newhistogram at 0, file u (f($2)) title col
plot newhistogram at 1, file u (f($3)) title col
plot newhistogram at 2, file u (f($4)) title col
set boxwidth 0.3 
plot newhistogram at 2.65, file  u (f($5)) title col
unset multiplot

This is somewhat better, for the colours are now consistent, and we also see that the last column has a different width. We also see how the positioning works: the right hand side of Japan's column is at 2.5, and since the width of Nauru's column is 0.3, its centre has got to be shifted by 0.15 with respect to 2.5. That adds up to 2.65. However, if we watch closely, we will also notice that the ytics and labels are drawn four times; after all, we have four plots. What, if we unset the ytics after the first plot? Well, we would end up with this

Rather upsetting! The problem is that once the tics are unset, the size of the figure changes, so we can no longer count on the plots' proper alignment. However, there is an easy remedy for this: all we have to do is not to unset the ytics, but to set them invisible. That is, we can do

plot newhistogram at 0, file u (f($2)) title col
set ytics ("      ", 30000)
plot newhistogram at 1, file u (f($3)) title col
plot newhistogram at 2, file u (f($4)) title col
set boxwidth 0.3 
plot newhistogram at 2.65, file  u (f($5)) title col

where we have 6 white spaces in the quote. You might wonder why on Earth 6. Well, the answer is that the label "30000" is actually " 30000", which takes up 6 characters' space. With this trick, we get

We have already achieved quite a lot, and slowly, but surely, we are getting to our goal. Just do not despair!

The next thing that we would need is proper scaling of the columns: we want all of them to be between 0 and 100 (%), i.e., we would have to sum all columns first, and then divide the values by the sum. And that is the snag: we have four columns, and we have to do the summing for each column independently, and before the final plots. Otherwise, our multiplot will be messed up. And this is where the array comes in handy: if we just had an array, and could retrieve values from it, we would be saved. And of course, we can do this. Let us take a small detour!

If we think about it, the array (5, 4, 6, 7, 8) is nothing but a finite series: its first element is 5, second element is 4, and so on. But we could also look at the series as a function: a mapping from the set of natural numbers to, well, to anything. In the example above, to natural numbers. It doesn't matter. My point is that an array is a function, a function for which h(0) = 5, h(1) = 4, h(2) = 6, h(3) = 7, and h(4) = 8. As long as this is true, we do not care what h(1.1) is. We need the function's values only at integer numbers. Then the only question is how we could define this function "on the fly". Being a physicist, and a lazy man, I would propose the following:

g(x,a) = (abs(x-a) < 0.1 ? 1 : 0)
h(x) = 5 * g(x,0) + 4 * g(x,1) + 6 * g(x,2) + 7 * g(x,3) + 8 * g(x,4)

g(x,a) is (apart from some numerical factors) nothing but a very primitive representation of a Dirac-delta, centred on 'a'. You can convince yourself that h(x) defined in this way fulfils the requirements above.

After this digression, let us see what we can do with this, and issue the following commands!

ARRAY = "h(x) = 0"
array(x, counter) = ( ARRAY.sprintf(" + %f*g(x,%d)", x/100.0, counter+1) )
ff(x, counter) = (($0 > 0 ? ARRAY = array(x, counter-1) : 1), total = total + x, x)
plot 'marimekko.dat' using 0:(ff($2, 2))

At this point, the variable ARRAY should look something like this

ARRAY = "h(x) = 0 + 91.630000*g(x,2) + 35.470000*g(x,2) + 77.220000*g(x,2) + 48.370000*g(x,2) + 34.410000*g(x,2)"

and if we evaluate it, the function value h(2) returns the sum of the numbers in the second column. (Apart from a factor of 100, of course.) Note that in order to take out the first line, which is the header, we have to use the condition

($0 > 0 ? ARRAY = array(x, counter-1) : 1)

which updates ARRAY only if we are processing the second record, at least. Also note that in order to get the sum of all columns, all we have to do is call this plot as many times as many columns there are. In the light of this, our next script could be this

reset

file = 'marimekko.dat'
col = 4

g(x,a) = (abs(x-a) < 0.1 ? 1 : 0)
ARRAY = "h(x) = 0"
array(x, counter) = ( ARRAY.sprintf(" + %f*g(x,%d)", x/100.0, counter) )
 
ff(x, counter) = (($0 > 0 ? ARRAY = array(x, counter) : 1), x)
plot for [i=2:col+1] 'marimekko.dat' using 0:(ff(column(i), i)) 

set xrange [-1:3]
set yrange [0:110]

eval(ARRAY);

set style data histograms
set style histogram columnstacked
set style fill solid border -1
set boxwidth 1.0

set multiplot
plot newhistogram at 0, file u ($2/h(2)) title col
set ytics ("    " 20)
plot newhistogram at 1, file u ($3/h(3)) title col
plot newhistogram at 2, file u ($4/h(4)) title col
set boxwidth 0.3
plot newhistogram at 2.65, file u ($5/h(5)) title col
unset multiplot

and this results in this figure

So, we are almost there: the columns are rescaled, and placed neatly next to each other. The only missing ingredient is the setting of the widths. But that is really easy: we only have to determine what the grand total is, and then scale the columns accordingly. Our script can, then, be modified as follows

reset
file = 'marimekko.dat'
col = 4

total = 0.0
g(x,a) = (abs(x-a) < 0.1 ? 1 : 0)
ARRAY = "h(x) = 0"
array(x, counter) = ( ARRAY.sprintf(" + %f*g(x,%d)", x/100.0, counter+1) )

ff(x, counter) = (($0 > 0 ? ARRAY = array(x, counter-1) : 1), total = total + x, x)
plot for [i=2:col+1] 'marimekko.dat' using 0:(ff(column(i), i)) 

set xrange [-0.3:1]
set yrange [0:110] 
eval(ARRAY);

set style data histograms
set style histogram columnstacked
set style fill solid border -1   

total = total / 100.0
position = 0.0

set multiplot
set boxwidth h(2)/total
plot newhistogram at position, file u ($2/h(2)) title col

set ytics ("    " 20)
set boxwidth h(3)/total; position = position + (h(2)+h(3))/total/2.0
plot newhistogram at position, file u ($3/h(3)) title col

set boxwidth h(4)/total; position = position + (h(3)+h(4))/total/2.0
plot newhistogram at position, file u ($4/h(4)) title col

set boxwidth h(5)/total; position = position + (h(4)+h(5))/total/2.0
plot newhistogram at position, file u ($5/h(5)) title col
unset multiplot

and this is what we wanted!

Adding labels to the rectangles is relatively easy: we could do the following

plot file using (position):(l($5)):5 with labels tc rgb "#ffffff"

where l(x) is a function that keeps track of the previous values of the column, and adds them as new values are processed. The definition of this function should be trivial.

The last thing that I would add here is that by using macros, we can tidy up the script: we no longer would need all those long and repetitive lines. In fact, we could also add another instruction to our 'ff' function, which would generate the plot command. The advantage of that is that in this way, we do not have to repeat the plot commands four times: we simply put that in our for loop, and then evaluate the resulting string. I discussed this trick in my last post, so, if you are interested in the details, you can look it up there.

Bending the arrows - "delaying" the plot

2010-04-26T12:01:00.000-07:00

The other day, I would have needed a couple of curved arrows on my plot, so I started to work out a method to get what I wanted. This, however, turned out to be rather interesting, so I thought that I would share the details with you.

First, we should just define what I mean by a curved arrow. Perhaps, the easiest way to define it is to show a plot, similar to this

In gnuplot, when one wants an arrow, one can invoke the following command:

set arrow from 0,0 to 1,1

or something similar. This will produce a straight arrow from (0,0) to (1,1). But what if we wanted to have an arrow, which is not straight. Well, in this case, we set a very short arrow, and draw a curve separately. The key to this is to set the arrow in such a way that it is tangential to the curve at the end point. It is easy to see that the following script would just do that

reset
unset key
eps = 0.001

set style arrow 1 head filled size screen 0.03, 15, 45 lt -1
cut(x,x1,x3) = ((x >= x1 && x <= x1 + (1.0-eps)*(x3-x1)) ? 1.0 : 1/0)
f(x) = 0.5+(x-1)*(x-1.2)*(x-1.4)

x1 = 0.5
x3 = 1.95
new_x = x1 + (1.0-eps)*(x3-x1)
set arrow from new_x, f(new_x) to x3,f(x3) as 1
plot [0:3] sin(x) with point ps 1 pt 6, f(x)*cut(x,x1,new_x) with line lt -1

First, we define an arrow style that we will use later. The arrow will be 0.03 screen sizes big, and the two angles determining the shape of the head are 15, and 45 degrees, respectively. Finally, we stipulate that the arrow be black, i.e. linetype -1. Then we define a window function, cut, which depends on the two end points, x1, and x3 (the reason for 3 will become clear soon), and the curve, f(x). In our plot, beyond what we actually want to plot, we will also plot f(x), but only between x1, and new_x, where new_x is a bit off with respect to the second end point. The degree of "bitness" is given by eps, which was defined at the beginning. However, before we actually plot anything, we have got to set the arrow, between new_x, f(new_x), and x3, f(x3). This construction ensures that the arrow is tangential to the curve.
At this point, we are ready to plot, which we actually execute in the next, and last line.

What we have created is great, but there are problems: first, we have to define our function, f(x), beforehand, we have to set the arrows by hand, and we also have to add the appropriate lines to our plot command. Quite tedious. There has got to be a better way!

For the say of example, let us suppose that we want a curved arrow that, say, connects (0,0) and (1,1) via a parabola that passes through the point (0.5, 0.25). If we are really pressed for it, we could do the following: First, we have to figure out the parameters of our parabola. In this case, it is quite easy, for it is nothing but x*x. Then we would draw a parabola between (0,0), and (0.99, 0.9801), and then draw an arrow from (0.99, 0.9801) to (1,1).

First, let us see, how we figure out the parameters of our parabola! We have two end points, and a "control" point, i.e., we have to solve the following set of equations

y1 = a*x1*x1 + b*x1 + c
y2 = a*x2*x2 + b*x2 + c
y3 = a*x3*x3 + b*x3 + c

for the unknown a, b, and c. You can convince yourself that the following will do

denom(x1, x2, x3) = x1*x1*(x2-x3) + x1*(x3*x3-x2*x2) + x2*x3*(x2-x3)
A(x1,y1,x2,y2,x3,y3) = ( (x2-x3)*y1 + (x3-x1)*y2 + (x1-x2)*y3 ) / denom(x1,x2,x3)
B(x1,y1,x2,y2,x3,y3) = ( (x3*x3-x2*x2)*y1 + (x1*x1-x3*x3)*y2 + (x2*x2-x1*x1)*y3 ) / denom(x1,x2,x3)
C(x1,y1,x2,y2,x3,y3) = ( (x2-x3)*x2*x3*y1 + (x3-x1)*x1*x3*y2 + (x1-x2)*x1*x2*y3 ) / denom(x1,x2,x3)
a = A(x1,y1,x2,y2,x3,y3)
b = B(x1,y1,x2,y2,x3,y3)
c = C(x1,y1,x2,y2,x3,y3)

We have done most of the hard work, the only thing that remains is how we "automate" this whole machinery, i.e., what do we do, if we have several arrows that we want to set. Again, as so many times in the past, we will utilise this new notion of function definition: the fact that a function is not only a x -> f(x) mapping, but this mapping, and a set of possibly unrelated instructions. What we will do is to define a "function" that sets our arrows, and, as the supplementary instruction, augments the plot command accordingly. First, let us take the following function definition

arrow(x1,y1,x2,y2,x3,y3) = (new_x = x1 + (1.0-eps)*(x3-x1), \
        a = A(x1,y1,x2,y2,x3,y3), b = B(x1,y1,x2,y2,x3,y3), c = C(x1,y1,x2,y2,x3,y3), \
        PLOT = PLOT.sprintf(", cut(x,%f,%f)*(%f*x*x+%f*x+%f) with lines lt -1", x1, x3, a, b, c), \
        ARROW.sprintf("set arrow from %f, %f to %f,%f as 1; ", new_x, a*new_x*new_x + b*new_x + c, x3, y3))

and try to understand what it does! For a start, it takes 6 arguments, which are nothing but the coordinates of the end points, and the control point. Then, it defines new_x, which we have already seen in the first example. In the next step, based on the 6 input arguments, calculates the three parameters of our parabola, and in the next line, adds the plot of this parabola to a string called PLOT. When adding to PLOT, we simply use the sprintf function. In the last line, we concatenate a string called ARROW, and another one, produced by another sprintf. It is easy to see that this sprintf returns the definition of an arrow between new_x, f(new_x), and x3, f(x3). We should also note that this line is the last line, which consequently means that whatever happens here is returned.

At this point we are really done, we only have to "populate" our plot. The full script takes on the form

reset
unset key
eps = 0.01

set style arrow 1 head filled size screen 0.03, 15, 45 lt -1

cut(x,x1,x3) = ((x >= x1 && x <= x1 + (1.0-eps)*(x3-x1)) ? 1.0 : 1/0)

denom(x1, x2, x3) = x1*x1*(x2-x3) + x1*(x3*x3-x2*x2) + x2*x3*(x2-x3)
A(x1,y1,x2,y2,x3,y3) = ( (x2-x3)*y1 + (x3-x1)*y2 + (x1-x2)*y3 ) / denom(x1,x2,x3)
B(x1,y1,x2,y2,x3,y3) = ( (x3*x3-x2*x2)*y1 + (x1*x1-x3*x3)*y2 + (x2*x2-x1*x1)*y3 ) / denom(x1,x2,x3)
C(x1,y1,x2,y2,x3,y3) = ( (x2-x3)*x2*x3*y1 + (x3-x1)*x1*x3*y2 + (x1-x2)*x1*x2*y3 ) / denom(x1,x2,x3)

ARROW = ""
PLOT = "p [0:3] sin(x) w p ps 1 pt 6"
arrow(x1,y1,x2,y2,x3,y3) = (new_x = x1 + (1.0-eps)*(x3-x1), \
        a = A(x1,y1,x2,y2,x3,y3), b = B(x1,y1,x2,y2,x3,y3), c = C(x1,y1,x2,y2,x3,y3), \
        PLOT = PLOT.sprintf(", cut(x,%f,%f)*(%f*x*x+%f*x+%f) with lines lt -1", x1, x3, a, b, c), \
        ARROW.sprintf("set arrow from %f, %f to %f,%f as 1; ", new_x, a*new_x*new_x + b*new_x + c, x3, y3))

ARROW = arrow(0,0,1,1.5,pi/2,1.03)
ARROW = arrow(0,0,1,0.3,pi/2,0.97)
eval(ARROW)
eval(PLOT)

which would result in the graph shown here:

Now it is clear what was PLOT: it is nothing, but the actual plot that we want to have. This is the string to which we concatenate our parabolae, one by one, every time we define a new arrow. After we defined all our arrows, we have two strings, ARROW, and PLOT. As such, they are no good, they will become instructions only when we evaluate them. That is what we do in the last two lines.

I would like to point out that my main reason for posting this was not that it can be used for creating curved arrows, but that this method is quite general. First, we can add to the plot, if that is needed, without having to keep track of all the tiny details. Second, the set command can be "fooled" by using the sprintf function. With the help of the string augmentation and the eval command, we can actually use parameters in our set instruction very efficiently.

Well, this is for today. I am waiting for suggestions as to what we should discuss next time. Cheers,
Zoltán

Bubble plots

2010-03-16T16:01:00.000-07:00

Yesterday, I discussed a method for adding an edge to an arbitrary symbol. If you recall (or roll down on this page), the idea was to trick gnuplot into plotting our data file twice, but in a way that each point was plotted twice in succession. Now, what if we plotted more times? There was really nothing special about the number 2, so there is no reason why we could not do this. But if we can, then we should, and see what comes out of it. With very small modifications, our script from yesterday can be turned into a bubble graph, like this

So, let us see how the machinery works!

reset
plot 'new_bubble1.dat' u 0:2
red_n = GPVAL_DATA_X_MAX

plot 'new_bubble2.dat' u 0:2
blue_n = GPVAL_DATA_X_MAX

plot 'new_bubble3.dat' u 0:2
green_n = GPVAL_DATA_X_MAX

rem(x,n) = x - n*(x/n)
size(x,n) = 3*(1-0.8*rem(x,n)/n)
c(x,n) = floor(240.0*rem(x,n)/n)
red(x,n) = sprintf("#%02X%02X%02X", 255, c(x,n), c(x,n))
blue(x,n) = sprintf("#%02X%02X%02X", c(x,n), c(x,n), 255)
green(x,n) = sprintf("#%02X%02X%02X", c(x,n), 255, c(x,n))

posx(X,x,n) = X + 0.03*rem(x,n)/n
posy(Y,x,n) = Y + 0.03*rem(x,n)/n

unset key
set border back
level = 40
plot for [n=0:level*(red_n+1)-1] 'new_bubble1.dat' using (posx($1,n,level)):(posy($2,n,level)) \
every ::(n/level)::(n/level) with p pt 7 ps size(n,level) lc rgb red(n,level) , \
for [n=0:level*(blue_n+1)-1] 'new_bubble2.dat' using (posx($1,n,level)):(posy($2,n,level)) \
every ::(n/level)::(n/level) with p pt 7 ps size(n,level) lc rgb blue(n,level) , \
for [n=0:level*(green_n+1)-1] 'new_bubble3.dat' using (posx($1,n,level)):(posy($2,n,level)) \
every ::(n/level)::(n/level) with p pt 7 ps size(n,level) lc rgb green(n,level)

Again, the first three plots are there for determining the sample size, and nothing more. We, thus, start out with a number of function definitions. The first one is a remainder function, the second one uses the remainder to return the size of the bubble, the third one is a simple helper function, returning values between 0 and 240, and red, blue, and green determine the colour of our bubbles. If you look carefully, you will notice that these colours are successively whiter as the remainder increases. Finally, again by making use of our remainder function, we define two position shifts: in order to give the impression that the bubbles are lit from the top right corner, we have to shift successive circles in that direction. The value of this shift is important in the sense that, if chosen too high, the circles belonging to the same data point will no longer cover each other. (This is not necessary a tragedy, see below.)

Then we decide to have 40 colour levels (we could have anything up to 255, although it might be a bit time consuming and unnecessary), and call our plots. The structure is the same as it was yesterday: we use a for loop for each data set, move the circles a bit, and set the colours to whiter shades. That is all.

Now, what happens, if we take too big a value for the shift? This, actually, might lead to interesting effects, as shown in this graph, where droplets represent the data points.

After having seen the simplest implementation, we should ask whether it is possible to add some decorations. E.g., whether it is possible to add a thin black edge to the symbols. It is relatively simple, as the following script shows. We only have to re-define some of our functions as follows

size(x,n) = (rem(x,n) == 0 ? 3.3 : 3*(1-0.8*rem(x,n)/n))
c(x,n) = floor(240.0*rem(x,n)/n)
red(x,n) = (rem(x,n) == 0 ? "#000000" : sprintf("#%02X%02X%02X", 255, c(x,n), c(x,n)))
blue(x,n) = (rem(x,n) == 0 ? "#000000" : sprintf("#%02X%02X%02X", c(x,n), c(x,n), 255))
green(x,n) = (rem(x,n) == 0 ? "#000000" : sprintf("#%02X%02X%02X", c(x,n), 255, c(x,n)))

posx(X,x,n) = (rem(x,n) < 2 ? X : X + 0.03*rem(x,n)/n)
posy(Y,x,n) = (rem(x,n) < 2 ? Y : Y + 0.03*rem(x,n)/n)

All these functions do is to check whether we are plotting the first round, and if so, set the colour to black. There is a small difference in the shifts, for we do not move the circles, if they are in the first or the second round. The reason is obvious, as is the result

OK, so we can plot bubbles, with or without black circumference, but we would also like to add a legend. Well, that is simple, in fact, nothing could be simpler. Just add the following the following three lines to our code

set label 1 'Red bubbles' at 9,6 left
set label 2 'Blue bubbles' at 9,5 left
set label 3 'Green bubbles' at 9,4 left

and the following six

for [n=0:level-1] 'new_bubble1.dat' using (posx(8.5,n,level)):(posy(6,n,level)) \
every ::(n/level)::(n/level) with p pt 7 ps size(n,level) lc rgb red(n,level) , \
for [n=0:level-1] 'new_bubble2.dat' using (posx(8.5,n,level)):(posy(5,n,level)) \
every ::(n/level)::(n/level) with p pt 7 ps size(n,level) lc rgb blue(n,level) , \
for [n=0:level-1] 'new_bubble3.dat' using (posx(8.5,n,level)):(posy(4,n,level)) \
every ::(n/level)::(n/level) with p pt 7 ps size(n,level) lc rgb green(n,level)

and we are done! All we do here is to plot our data files in a silly way: we plot a single point at (8.5,6), (8.5,5), and (8.5,4). The plotting of the data file does not happen in this sense, we use it for convenience's sake only. (This trick can also be used for the post from yesterday.) There, you have it!

Defining new symbols

2010-03-16T14:02:00.000-07:00

Some time ago, I showed a method with which we could add a "frame" to a symbol. If you recall, what we did was to plot everything twice, and in order to duplicate our data set, we used a simple gawk script. Now, there is another way of doing this, one which does not rely on the gawk script, in fact, on any external script. I will discuss this method today. The gist of the trick is discussed in the old post, therefore, you are encouraged to cast, at least, a cursory glance at that, if you haven't yet done it.

As I have already pointed out, we had to duplicate our data set. To be more accurate, we haven't got to duplicate anything, we have simply got to plot the data twice. Now, the difficulty is that is we do this in a primitive way, issuing the plot command twice, and taking the same data set, the points might overlap, and leads to some undesired results. So, the task is to plot the data set twice, but to plot each plot twice, and not the data set as a whole. For this, we will use the for loop introduced in gnuplot 4.4, and the 'every' keyword. To cut a long story short, I give my script here, and discuss it afterwards.

reset 
plot 'new_symbol1.dat' u 0:2
red_n = GPVAL_DATA_X_MAX

plot 'new_symbol2.dat' u 0:2
blue_n = GPVAL_DATA_X_MAX

plot 'new_symbol3.dat' u 0:2
green_n = GPVAL_DATA_X_MAX

parity(n) = (n/2.0 == int(n/2.0) ? 0 : 1)
size(n) = 2 - parity(n)*0.4
colour(n,r,g,b) = sprintf("#%02X%02X%02X", parity(n)*r, parity(n)*g, parity(n)*b)

unset key
set border back
plot for [n=0:2*red_n+1] 'new_symbol1.dat' using 1:2 \
every ::(n/2)::(n/2) with p pt 7 ps size(n) lc rgb colour(n,255,0,0) ,\
for [n=0:2*blue_n+1] 'new_symbol2.dat' using 1:2 \
every ::(n/2)::(n/2) with p pt 9 ps size(n) lc rgb colour(n,100,100,255) ,\
for [n=0:2*green_n+1] 'new_symbol3.dat' using 1:2 \
every ::(n/2)::(n/2) with p pt 5 ps size(n) lc rgb colour(n,0,150,0)

Then, let us see what we have here! The first 6 lines are only to retrieve the number of data points in our data sets. If you know this from somewhere else, you can skip these, with the caveat that 'red_n', 'blue_n', and 'green_n' should still be defined somewhere.

Next we define three functions, the first of which determines the parity of an integer, returning 1, if the number is odd, and 0, if it is even. The second function returns a number, depending on the parity of its argument. Surprising as it is, this function will determine the size if the symbol, when we plot. Finally, the third function returns a string, which is equal to the colour given by the triplet (r,g,b), if the first argument, 'n', is odd, and black, if the first argument is even. At this point, it should be clear that we could have defined a function that returns a different colour for even numbers.

We are done with everything, but the plotting, so let us do that! As you see, for each data set, we step through the numbers, but not once, but twice: first plotting in black, and second, plotting with some decent colour. At the same time, we change the symbol size, so that the black symbols are always a bit bigger, than the red, blue, or green. Once all three plots have been called, the following graph will appear:

We can see that the symbols overlap each others, as they should. Now, what about the keys, should we need them? Well, that requires some handwork, but it is not hard, actually. The following self-explanatory script should do

set label 1 'Red symbols' at 1.3, 8 left
plot for [n=0:2*red_n+1] 'new_symbol1.dat' using 1:2 \
every ::(n/2)::(n/2) with p pt 7 ps size(n) lc rgb colour(n,255,0,0), \
n=0, '-' using 1:2 with p pt 7 ps size(n) lc rgb colour(n,255,0,0), \
n=1, '-' using 1:2 with p pt 7 ps size(n) lc rgb colour(n,255,0,0)
1 8
e
1 8
e

and this produces the following figure

Phong on histograms with a one-liner (almost:-)

2010-02-26T13:13:00.000-08:00

We have seen in the last couple of posts that with the new concept of functions, quite a few interesting effects can be achieved. Today I would like to show a trick that solves a problem that I discussed some time ago, when we made shiny histograms using a for loop in gnuplot. We will do the same thing here, but in two lines only. It is quick, and the results are just as good as in that case.

So, here is my data file

which I will just name as 'bar.dat', and here is our script

reset
unset key
set style fill solid 1.0
set yrange [0:6]

colour = "#080000"
f(x,n) = (colour = sprintf("#%02X%02X%02X", 128+n/2, n, n), x)
w(n) = 0.8*cos(n/230.0*pi/2.0)

plot for [n=1:230:2] 'bar.dat' u 0:(f($1,n)):(w(n)) with boxes lc rgbcolor colour

Simple enough, let us see what it does! The first four lines are just the usual settings, although, the yrange is really irrelevant. I set it only for aesthetic reasons (otherwise, gnuplot would set the yrange automatically to [1:5] for the data file above, and we wouldn't see one of the columns). Then we define a variable called 'colour'm which we will eventually overwrite in our function definition of f(x,n). f(x,n) returns x, thus, in this regard it would be absolutely useless, but when doing so, it actually prints a string to 'colour'. The next function is w(n), which will determine in what fashion our colour will converge to white.
Finally, we plot the data file some 115 times, each time with a smaller, and shinier box. At the end, we get something like this

We can very easily change the direction of the light. All we have to do is define a new function that shifts the bars as we progress with our for loop. So, the new script could be something like this

reset
unset key
set style fill solid 1.0
set yrange [0:6]

colour = "#080000"
f(x,n) = (colour = sprintf("#%02X%02X%02X", 128+n/2, n, n), x)
w(n) = 0.8*cos(n/230.0*pi/2.0)
shift(x,n) = x-0.8*n/850.0

plot for [n=1:230:2] 'bar.dat' u (shift($0,n)):(f($1,n)):(w(n)) with boxes lc rgbcolor colour

with a result as in this graph

It should be really easy to modify the script to accommodate more data sets. Well, this is for now. I don't actually know what I will write about next time, but I am sure that there will be something!
Cheers,
Zoltán

Parametric plot from a file II.

2010-02-25T13:38:00.000-08:00

As I promised yesterday, we will take a closer look at the pie chart, once more, and see how we can utilise what we have learnt recently. I should point out here, that this is not the only way of plotting a pie from a file. If you feel like building your gnuplot from source, you can check out either the CVS tree, or the patch tracker, where you can find a patch that makes it possible to plot slices. You can see a demo here. But we will try a different route here.

First, here is our data file (it could be anything, really)

1 1 Dolphins
2 1 Whales
2 0 Sharks
3 0 Penguins
4 1 Kiwis
5 0 Tux

and here is our script

reset
unset key; set border 0; unset tics; unset colorbox; set size 0.6,1
set urange [0:1]
set vrange [0:2*pi]
set macro
sum = 0.0
ssum = 0.0
n = 0
PLOT = "splot 0, 0, 1/0 with pm3d"

count(x) = (ssum = ssum + $1, 1)
g(x,y,n) = \
sprintf(", \
u*cos((%.2f+%.2f*v)/ssum), \
u*sin((%.2f+%.2f*v)/ssum), \
%d @PL", x, y, x, y, n)

f(x) = (PLOT = PLOT.g(2*pi*sum, x, n), sum = sum+x, n = n + 1, x)

plot 'new_pie.dat' u 1:(count($1)), '' u 1:(f($1))

PL = "with pm3d"
set parametric; set pm3d map; 
eval(PLOT)

There is really nothing that we haven't discussed before: we set a couple of things at the beginning, but most importantly, the macro, and sum, ssum, and n. Then we define a string, PLOT, and two functions. One is to sum the values in our file (we need this, so that we can scale the full range of angles to two pi), and another one, that writes our PLOT command for later use. Note that the first plot in PLOT is actually empty, we plot 1/0. This seems a bit silly, doesn't it? Well, it does, but there is a good reason: successive plots must be separated by a comma, and if we have an empty plot at the very beginning, then we can put the commas before the plots, not after, and in this way, we needn't keep track of which plot we are actually processing. Remember, yesterday we used a separate counter, and an if statement, to determine, whether we need the comma, or not. This we can avoid here.
Next we call the two dummy plots, and finally, we evaluate our PLOT string. Oh, no! At the very end, we marvel in awe at the figure that we produced.

So far, so good, but what if we wanted to add labels, e.g., the value of the slice? That is really easy. All we have to do is to define a function that produces the label. Here is our updated script

reset
unset key; set border 0; unset tics; unset colorbox; set size 0.6,1
set urange [0:1]
set vrange [0:2*pi]
set macro
sum = 0.0
ssum = 0.0
n = 0
PLOT = "splot 0, 0, 1/0 with pm3d"
LABEL = ""

count(x) = (ssum = ssum + $1, 1)
g(x,y,n) = \
sprintf(", \
u*cos((%.2f+%.2f*v)/ssum), \
u*sin((%.2f+%.2f*v)/ssum), \
%d @PL", x, y, x, y, n)

lab(alpha, x) = sprintf("set label \"%s\" at %.2f, %.2f; ", \
x, 1.2*cos(alpha), 1.2*sin(alpha))

f(x) = (PLOT = PLOT.g(2*pi*sum, x, n), \
LABEL = LABEL.lab(2*pi*sum/ssum+pi*x/ssum, sprintf("%2.f", x)), \
sum = sum+x, n = n + 1, x)

plot 'new_pie.dat' u 1:(count($1)), '' u 1:(f($1))

PL = "with pm3d"
set parametric; set pm3d map; set border 0; unset tics; unset colorbox;
set size 0.6,1
eval(LABEL)
eval(PLOT)

We have an addition string, LABEL, which we initialise with the value "". Then we define a function that prints "set label ..." with the proper positions, and finally, we insert this function in the definition of f(x). Of course, once we called f(x) in the plot, we have to evaluate the string LABEL. So, this is what we get

This script can trivially be modified to print strings that are stored in our file. Watch just the following two lines

...
lab(alpha, x) = sprintf("set label \"%s\" at %.2f, %.2f centre; ", \
x, 1.2*cos(alpha), 1.2*sin(alpha))

f(x) = (PLOT = PLOT.g(2*pi*sum, x, n), \
LABEL = LABEL.lab(2*pi*sum/ssum+pi*x/ssum, stringcolumn($3)), \
sum = sum+x, n = n + 1, x)
...

and we are done. Here is the new pie

Now, you might wonder why we had three columns in our data file, if we didn't want to use it. Well, perhaps, we wanted, just haven't got time till now. So, what could we do with those ones and zeros in the second column? We will make the pie explode! It is really simple, we have to modify two lines in our last script

...
g(x,y,n,dx,dy) = \
sprintf(", \
%.2f+u*cos((%.2f+%.2f*v)/ssum), \
%.2f+u*sin((%.2f+%.2f*v)/ssum), \
%d @PL", 0.2*dx, x, y, 0.2*dy, x, y, n)

lab(alpha, x, r) = sprintf("set label \"%s\" at %.2f, %.2f centre; ", \
x, (1.25+0.2*r)*cos(alpha), (1.25+0.2*r)*sin(alpha))

f(x) = (PLOT = PLOT.g(2*pi*sum, x, n, \
$2*cos(2*pi*sum/ssum+pi*x/ssum), \
$2*sin(2*pi*sum/ssum+pi*x/ssum)), \
LABEL = LABEL.lab(2*pi*sum/ssum+pi*x/ssum, stringcolumn($3), $2), \
sum = sum+x, n = n + 1, x)
...

And here is the pie, when exploded

I should mention here that if you are not happy with the colours, it is really easy to help it: all we have to do is to modify the colour palette, using whatever colour combinations. We have covered a lot of material today. Till next time,
Gnuplotter

Plotting in 6 dimensions - parametric plot from a file

2010-02-24T15:01:00.000-08:00

Today I would like to touch on a vast subject, so prepare for a long post. However, I hope that the post will be worthwhile, for I want to discuss something that cannot be done in any other way. In due course, we will see how we can use gnuplot to create parametric plots from a file. What I mean by that is the following: if you want to plot, say, 10 similar objects, whose size is determined by the first column in a file. Of course, there are cases, when one can manipulate the size, e.g., if there is a pre-defined symbol, we can use one of the columns in a file to determine the size of the symbol. As an example, we can do this

plot 'foo' u 1:2:3 with point pt 6 ps var

which will draw circles whose radius is given by the third column in 'foo'. This is all well, but it works for a limited number of cases only, namely, when there is a symbol to start out with. But what happens, if we want to draw an object that is not a symbol, e.g., arcs of a circle, whose angle is given by one of the columns in a file, or cylinders, whose height is a variable, read from a file. As you can guess from these two suggestions, what we will do is to draw a pie, and a bar chart. I understand that we have done this a couple of times before, but this time, we will stay entirely in the realm of gnuplot, and the scripts are really short. We just have to figure out what to write in the scripts. But beyond this, I will also show how we can plot in 6 dimensions. We will plot ellipses on a plane (first 2 columns), whose two axes are given by the 3rd and 4th axis, the orientation by the 5th, and the colour by the 6th. If you are really pressed for it, you can add three more dimensions: if you draw ellipsoids in 3D, take all three axes from a file, and also the orientation, that would make 9 dimensions altogether. Quite a lot!

So, let us get down to business! The first thing that I would like to discuss is the evaluate command. This is a really nifty way of shortening repetitive commands. Let us suppose that we want to place 10 arrows on our graph, and only the first coordinate of the arrows changes, otherwise everything is the same. Setting one arrow would read as follows

set arrow from 0, 0 to 1, 1

Of course, there are quite a few settings that we could specify, but this was supposed to be a minimal example. Then, the next arrow should be

set arrow from 1, 0 to 2, 1

and so on. What if we do not want to write this line a thousand times, and we do not want to search for the coordinate that we are to change, the first one, in this case? We could try the following

a(x) = sprintf("set arrow from %d, 0 to %d, 1", x, x+1)

This function takes 'x', and returns a string with all the settings and coordinates. So, we are almost done. The only thing we should do is to make gnuplot understand that what we want it to treat a(x) as a command, not as a string. Enter the eval command: it takes whatever string is presented to it, and turns it into a command. Thus, the following script creates 5 arrows, all parallel to each other, and consecutively shifted to the rigth

a(x) = sprintf("set arrow from %d, 0 to %d, 1", x, x+1)
eval a(0)
eval a(1)
eval a(2)
eval a(3)
eval a(4)

I believe, this is a much simpler and cleaner procedure, than this

set arrow from 0, 0 to 1, 1
set arrow from 1, 0 to 2, 1
set arrow from 2, 0 to 3, 1
set arrow from 3, 0 to 4, 1
set arrow from 4, 0 to 5, 1

I should mention here that if chunks of a command are the same, another method of abbreviating them is to use macros. Those are disabled by default, so first we have to set it. Then it works as follows

set macro
ST = "using 1:2 with lines lt 3 lw 3"
plot 'foo' @ST, 'bar' @ST

i.e., the term @ST is expanded using the definition above, therefore, this plot is equivalent to this one

plot 'foo' using 1:2 with lines lt 3 lw 3, 'bar' using 1:2 with lines lt 3 lw 3

but the previous one is much more readable. I would also say that using capitals for the macros is probably not a bad idea, because then they cannot be mistaken for standard gnuplot commands. This much in the way of macros!

So, we have the evaluate command, and we have a new concept for functions. Then let us take a closer look at the following code

a(x) = sprintf("set arrow from %d, 0 to %d, 1;\n", x, x+1)
ARROW = ""
f(x) = (ARROW = ARROW.a(x), x)
plot 'foo' using 1:(f($1))

and let us suppose that our file 'foo' contains the following 5 lines

After plotting 'foo', the string 'ARROW' will be the following

set arrow from 1, 0 to 2, 1;
set arrow from 3, 0 to 4, 1;
set arrow from 5, 0 to 6, 1;
set arrow from 7, 0 to 8, 1;
set arrow from 9, 0 to 10, 1;

I.e., we have a string, which contains instructions for setting 5 arrow. If, at this point, we simply evaluate this string, all 5 arrows will be set. Therefore, we have found a way of using a file to set the coordinates of an arrow. (N.B., if it was for the arrows only, we wouldn't have had to do anything, since there is a plotting style, 'with vector', as we discussed some weeks ago.)

We will use this trick to create a parametric plot, taking parameter values from a file, first plotting the ellipses! Again, we have got to create some dummy data, and since we now need 6 columns, we will use the errorbars

reset
f(x) = rand(0)
set sample 50
set table 'ellipse.dat'
plot [0:10] '+' using (20*f($1)):(20*f($1)):(f($1)):(f($1)):(3.14*f($1)):(f($1)) w xyerror
unset table

which will produce 6 columns and 50 lines. Having produced some data, let us see what we can do with it. Here is our script:

PRINT(x, y, a, b, alpha, colour) = \
sprintf("%f+v*(%f*cos(u)*cos(%f)-%f*sin(u)*sin(%f)),
%f+v*(%f*cos(u)*sin(%f)+%f*sin(u)*cos(%f)),
%f with pm3d", x, a, alpha, b, alpha, y, b, alpha, b, alpha, colour)
PLOT = "splot "
num = -1
count(x) = (num = num+1, 1)
g(x) = (PLOT = PLOT.PRINT($1, $2, $3, $4, $5, $6), \
($0 < num ? PLOT=PLOT.sprintf(",\n") : 1/0))
plot 'ellipse.dat' u 1:(count($1))
plot 'ellipse.dat' using 1:(g($1))

unset key
set parametric
set urange [0:2*pi]
set vrange [0:1]
set pm3d map
set size 0.5, 1
eval(PLOT)

First, we have the definition of a print function that looks rather ugly, but is quite simple. We want to plot

a*v*cos(u), b*v*sin(u), colour

where a, and b are the axes of the ellipse, and colour is going to specify, well, its colour. However, we want to translate the ellipse to its proper position, and we also want to rotate it by an amount given by the 5th column, so we have to apply a two-dimensional rotation on the object. Therefore, we would end up with a function similar to this

x+v*(a*cos(u)*cos(alpha)-b*sin(u)*sin(alpha)), y + v*(a*cos(u)*sin(alpha)+b*sin(u)*cos(alpha)), colour

Now you know why that print function looked so complicated! After this, we define a string, PLOT, that we will expand as we read the file. But before that, we have to count the lines in the file. The reason for that is that successive plots must be separated by a comma, but there shouldn't be a comma after the last plot. So, we just have to know where to stop placing commas in our string. Then we define the function that does nothing useful, but concatenates the PLOT string as it reads the file. Here we use the number of lines that we determine in a dummy plot. At this point we are done with the functions, all we have to do is plotting.
First we count, then plot g(x). At this point, we have the string that we need. We only have to set up our plot. Remember, we have a parametric plot, where the range of one of the variables is in [0:2*pi], while the other one is in [0:1]. Easy. Then we just have to evaluate our plot string, and we are done. Look what we have made here: a six-dimensional plot!

I think that this script is much less complicated, than many that we have discussed in the past. Short and clear, thanks to the eval command, and the new concept of functions. Besides, we pulled off a trick that was impossible by other means. I started out saying that we will create bars and pie. I believe, having seen the trick, it should be quite simple now, but in case you insist on seeing it, I will discuss it in my next post.

The map, the inline function, and the macro

2010-02-10T13:23:00.000-08:00

Some time ago, on the 26th of July, to be more accurate, I showed how somewhat decent-looking maps can be created with gnuplot. With the wisdom of hindsight, that was a rather ugly hack, I must say. Even worse, it seems that it is not quite fail-safe. At least, I have obtained reports complaining about it. Could we, then, do better? Could we, perhaps, throw out that disgusting gawk script, with all the hassle that comes with it? Could we, possibly, manage the whole affair in gnuplot? Sure, we could. And here is how, just keep reading!

On the 17th of January, we saw that in the new version of gnuplot, functions take on a funny property, namely, they can contain algebraic statements not related to the return value. We also saw that this feature can be used to perform searches of some sort: as we "plot" a file, and step through the numbers in the file, we can assign values to variables, provided that some conditions are fulfilled. It is easy to see that in this way, we can determine the minimum or the maximum of a data file, e.g. But we can do much more, than that.

We should also recall from that old post on the map what the contour file looks like. In case you have forgotten, here is a small section of it

# Contour 0, label:        2
-0.391812  3.63636  2
...
-0.959596  3.50978  2

-0.959596  3.50978  2
...
-0.391812  3.63636  2


# Contour 1, label:      1.5
-1.20098  4.51515  1.5
-1.16162  4.54423  1.5
-1.15982  4.54545  1.5
...

What we have to realise is the following: first, contours lines belonging to the same level are not necessarily contiguous (this is quite obvious, for there is no reason why they should be), and if there is a discontinuity, it manifests itself in a single blank line in the contour file, and second, contour lines belonging to different levels are separated by two blank lines. So, in the data file above, there is a blank between the lines -0.959596 3.50978 2, and
-0.959596 3.50978 2, and there are two blanks between -0.391812 3.63636 2, and # Contour 1, label: 1.5. By the way, the third column is the value of that particular contour line.

This observation has at least one important consequence: we can decide which contour line we want to plot, simply by using the index keyword. You might recall, that indexing the data file pulls out one data block, which is defined by a chunk of data flanked by two blank lines.

Now, what about the labels, and the white space that they need? Well, the white space is quite easy: what we will plot is not the contour line, but a function, which returns an undefined value at the place of the white space, e.g., this one (whatever eps and xtoy mean)

f(x,y) = ((x-x0)*(x-x0)+(y-y0)*(y-y0)*xtoy*xtoy > eps ? y : 1/0)

Normally we would plot the contour lines as

plot 'contour.dat' using 1:2 with lines

but instead of this, now we will use this

plot 'contour.dat' using 1:(f($1,$2)) with lines

This will leave out those points which are too close to (x0, y0). And the labels? Well, that is not difficult either. Take this function

lab(x,y) = ( (x == x0 && y == y0) ? stringcolumn(3) : "")

and this plot

plot 'contour.dat' using 1:2:(lab($1,$2)) with labels

This will put the labels at (x0, y0), and even better, we haven't got to set the labels by hand, they are taken from the data file.

So, we have seen how we can plot the contour, leave out some white space, and then put a label at that position. The only remaining question is how we determine where the label should be. And this is where we come back to our inline functions. For the sake of example, let us take this function and the accompanying plot

g(x,y)=(((x > xl && x < xh && y > yl && y < yh) ? (x0 = x, y0 = y) : 1), 1/0)
plot 'contour.dat' using 1:(g($1,$2))

What on Earth does this plot do? The plot itself does absolutely nothing: it is always 1/0. However, while we are doing this, we set the value of x0, and y0, if the two arguments are not too close to the edge of the plot. This latter condition is needed, otherwise labels could fall on the border, which doesn't look particularly nice.

By now, we have all the bits and pieces, we have only got to put them together. Let us get down to business, then!

I will split the script into two: the first produces the dummy data, while the second does the actual plotting. So, first, the data production.

reset
filename = "cont.dat"
xi = -5; xa = 0; yi = 2; ya = 5;
xl = xi + 0.1*(xa - xi); xh = xa - 0.1*(xa-xi);
yl = yi + 0.1*(ya - yi); yh = ya - 0.1*(ya-yi);
xtoy = (xa-xi) / (ya-yi)
set xrange [xi:xa]
set yrange [yi:ya]
set isosample 100, 100
set table 'test.dat'
splot sin(1.3*x)*cos(.9*y)+cos(.8*x)*sin(1.9*y)+cos(y*.2*x)
unset table
set cont base
set cntrparam level incremental -3, 0.5, 3
unset surf
set table filename
splot sin(1.3*x)*cos(0.9*y)+cos(.8*x)*sin(1.9*y)+cos(y*.2*x)
unset table

What we should pay attention to here is the definition of a handful of variables at the very beginning. Some are already obvious, like xi, xa and the like, and some will become clear in the second part. Now, the plotting takes place here

reset
unset key
set macro
set xrange [xi:xa]
set yrange [yi:ya]

set tics out nomirror
set palette rgbformulae 33,13,10
eps = 0.05

g(x,y)=(((x > xl && x < xh && y > yl && y < yh) ? (x0 = x, y0 = y) : 1), 1/0)
f(x,y) = ((x-x0)*(x-x0)+(y-y0)*(y-y0)*xtoy*xtoy > eps ? y : 1/0)
lab(x,y) = ( (x == x0 && y == y0) ? stringcolumn(3) : "")

ZERO = "x0 = xi - (xa-xi), y0 = yi - (ya-yi), b = b+1"
SEARCH = "filename index b using 1:(g($1,$2))"
PLOT = "filename index b using 1:(f($1,$2)) with lines lt -1 lw 1"
LABEL = "filename index b using 1:2:(lab($1,$2)) with labels"

b = 0
plot 'test.dat' with image, \
@SEARCH, @PLOT, @LABEL, @ZERO, \
@SEARCH, @PLOT, @LABEL, @ZERO, \
@SEARCH, @PLOT, @LABEL, @ZERO, \
@SEARCH, @PLOT, @LABEL, @ZERO, \
@SEARCH, @PLOT, @LABEL, @ZERO, \
@SEARCH, @PLOT, @LABEL, @ZERO, \
@SEARCH, @PLOT, @LABEL, @ZERO, \
@SEARCH, @PLOT, @LABEL, @ZERO, \
@SEARCH, @PLOT, @LABEL, @ZERO, \
@SEARCH, @PLOT, @LABEL, @ZERO, \
@SEARCH, @PLOT, @LABEL, @ZERO, \
@SEARCH, @PLOT, @LABEL

A bit convoluted, isn't it? OK, we will walk through the script, line by line.

First is the range setting, and then ticks go to the outside, just for aesthetic reasons. We also define eps here, which determines how much "white space" we have for the labels. Then we define the three functions that we discussed above. We have already seen eps, and the meaning of it, but what about xtoy? Despite its name, this is not something to play with, rather the ratio of x to y, or more precisely, the ratio of the xrange to the yrange. This is needed, if the two ranges are of different order of magnitude, e.g., if xrange is something like [0:1], while yrange is [0:1000]. But this ratio is automatically calculated at the beginning, you haven't got to worry about it.

After this, we define 4 macros. These are abbreviations for longer chunks of code, and make life really easier. The idea is that when confronted with a macro, gnuplot expands it as a string, and then acts accordingly. In my opinion, if written properly, macros can make the script rather readable.

The first of the macros, ZERO, is needed, because in our SEARCH macro, which is nothing but a call to the function g(x,y), if the condition is not satisfied for a particular data block, then x0, y0 wouldn't be updated, therefore, the label would end up at the wrong position. At the same time, ZERO also increments the value of b, which determines which data block we are actually plotting. b is used in the indexing in the macros SEARCH, PLOT, and LABEL. We have already mentioned SEARCH, PLOT plots the contour with the white space at the position given by x0, and y0 (this is calculated in the SEARCH macro), and finally, LABEL places the value of the contour line at that position.

At this point, we have defined everything, all that is left is plotting. We do it 13 times, because our zrange, or the contour lines were given between -3, and 3, with steps of 0.5. In this particular case, there are only 10 contour lines, and gnuplot will complain that the last 3 data blocks are empty, but this is not an error, only a warning. Shouldn't we look at the figure, perhaps? But of course! Here it is:

The only thing that I should like to point out is that the white space is made for a particular contour line, but there is no guarantee that, if the contour lines are too close to each other, the label does not cover a neighbouring contour line. If that happens, I would simply suggest to increase the contour spacing by incrementing the parameter in the set cntrparam line.

I hope that this method proves better, than the other one, and that it will be easier to use. In the next post, I will re-visit the inline functions, and show a nifty trick with them. Cheers,
Gnuplotter

Further new features in gnuplot 4.4

2010-01-17T13:38:00.000-08:00

Today I would like to discuss another new feature in gnuplot 4.4. This will be the notion of "inline" data manipulation. I don't really know what the proper name would be for this feature, so I will just show what it is.

Normally, when one plots a function or file, the command has the following structure

plot 'foo' using 1:2 with line, f(x) with line

That was the old syntax. In the new version of gnuplot, we can insert arithmetic expressions in the plot command as follows

f(x) = a*sin(x)
plot a = 1.0, f(x), a = 2.0, f(x)

Now, this has some implications. First, one has to be a bit careful, because the arithmetic expressions are separated from the actual function by a comma. However, the 'for' loop that we discussed a week ago, reads statements up to the comma, and then returns to the beginning of the statement. In other words,

plot for [i=1:10] a = i, f(x)

will evaluate the expression a = i ten times, and then plots f(x). At that point, the value of 'a' will be 10, therefore, we have only one plot, and that will be 10*sin(x).

The second implication is that the notion of a function has completely changed. What we do in a plot command now is no longer a mapping of the form

x -> f(x)

but rather, the evaluation of a set of instructions, one of which is the above-mentioned mapping. But the crucial point here is that the mapping is not the only allowed statement. The upshot is that "functions" have become a set of operations, and the following statement is completely legal

f(x) = (a = a+1.0, a*sin(x))
a = 0.0
plot for [i=1:10] f(x)

(It is a complete different question whether this plot makes any sense...)

What we should notice is the fact that now a function can have the form of

f(x) = (statement1, statement2, statement3, return value)

and when the function is called, statement1, statement2, statement3 are evaluated, and 'return value' is returned. We should not underestimate the significance of this! Many things can be done with this. I will show a few of them below.

The first thing I would like to dive into is calculating some statistics of a file. Let us see how this works!

reset
set table 'inline.dat'
plot sin(x)
unset table
num = 0
sum = 0.0
sumsq = 0.0

f(x) = (num = num+1, sum = sum+x, sumsq = x*x, x)
plot 'inline.dat' using 1:(f($2))

print num, sum, sumsq

which will print out

100 -1.77635683940025e-15 0.295958848441

We expected this, for the number of samples is 100 by default, and the sum should be 0 in this case.

So, what about finding the minimum and its position in a data file? This is quite easy. All we have to do is to modify our function definition, and insert a statement that determines whether a value is minimal or not.

reset
set table 'inline.dat'
plot sin(x)
unset table

num = 0
min = 1000.0
min_pos = 0
min_pos_x = 0

f(x,y) = ((min > y ? (min = y, min_pos_x = x, min_pos = num) : 1), num = num+1,  y)
plot 'inline.dat' using 1:(f($1, $2))

print num, min, min_pos_x, min_pos

which prints

100 -0.999385 4.74747 73

i.e., the minimum is at the 73rd record (we count from 0), at x = 4.74747, and its value is -0.999385. Note that instead of an 'if' statement, we use the ternary operator to decide whether min, min_pos_x, and min_pos should be updated.

The implementation of calculating the standard deviation, e.g., should be trivial:

sum = 0.0
sumsq = 0.0
f(x) = (num = num+1, sum = sum + x, sumsq = sumsq + x*x, x)
plot 'inline.dat' using 1:(f($1))

print num, sqrt(sumsq/num - (sum/num)*(sum/num))

We have thus seen how the "inline" arithmetic can be used for calculating quantities, e.g., various moments, minima/maxima and their respective positions. These involve the sequential summing or inspection of the data set. But this trick with the function definition can be used for back-referencing, too. This is what we will discuss next.

The trick is to use a construct similar to this

backshift(x) = (prev = pres, pres = x, prev)

which will store the last but one value in the variable 'prev', and return it. That is, the following code shift the whole curve to the right by one

reset
set table 'inline.dat'
plot sin(x)
unset table

pres = 0.0

backshift(x) = (prev = pres, pres = x, ($0 > 0 ? prev : pres))
plot 'inline.dat' using 1:(backshift($2)) with line, '' u 1:2 with line

(In cases like this, we always have to decide what to do with the first/last data record. In this particular case, I opted for duplicating the first record, - this is what happens in the ternary operator - but this is not the only possibility.) If, for some reason, you have to shift the curve by more, you do the same thing, but multiple times. E.g., the following code shifts by 3 places.

backshift(x) = (prev1 = prev2, prev2 = prev3, prev3 = x, prev1)

Once we have this option of back-referencing, we should ask the question what it can be used for. I show two examples for this.
The first example is drawing arrows along a line given by the data set. Drawing arrows one by one is done by using

set arrow from x1,y1 to x2,y2

but we have to use a different method, if we want to plot the arrows from a file. Incidentally, there is a plotting style, 'with vectors', that works as

plot 'foo' using 1:2:3:4 with vectors

where the first two columns specify the coordinates of the beginning, and the second two columns specify the relative coordinates of the vectors. So, it works on four columns. What should we do, if we want to plot vectors from the points in a file. Well, we use the back shift that we defined above. Our script is as follows:

reset
unset key
set sample 30
set table 'arrowplot.dat'
plot [0:3] sin(x)+0.2*rand(0)
unset table

px = NaN
py = NaN
dx(x) = (xd = x-px, px = ($0 > 0 ? x : 1/0), xd)
dy(y) = (yd = y-py, py = ($0 > 0 ? y : 1/0), yd)

plot 'arrowplot.dat' using (px):(py):(dx($1)):(dy($2)) with vector

which results in the following figure:

Note that we used the ternary operator to get rid of the very first data point. This is needed, because the arrows connect two points, that is, there will be one less arrow, than data points.

In the second example, we will turn this around. In my post in last August, plotting the recession, I showed how the background of a plot can be changed, based on whether the the curve is dropping, or increasing. Let us take the following script

reset
set sample 20
set table 'inline.dat'
plot [0:10] exp(-x)+1.0+rand(0)
unset table

unset key

px = 0
py = 1000
dx(x) = (xd = x-px, px = x, xd)
dyn(y) = (yd = y-py, py = y, (yd < 0 ? yd : 1/0))
dyp(y) = (yd = y-py, py = y, (yd >= 0 ? yd : 1/0))

plot 'inline.dat' using (px):(py):(dx($1)):(dyp($2)) with vector nohead lt 1 lw 3, \
px = 0, py = 0, '' using (px):(py):(dx($1)):(dyn($2)) with vector nohead lt 3 lw 3

which creates the following graph

First we produce some data; old trick. Then we take our difference functions, in this three of them. The first one is identical to that in the previous script. The second and the third are identical, except that the second returns a sensible value, if and only, if the slope is negative, while the third one returns 1/0, if the slope is negative. Then we just plot our data, making sure that we re-initialise px, and py before the second plot. Simple.

Another utilisation of the back reference can be found on gnuplot's web site, under running averages.

Next time I will try to go a bit further, and demonstrate some other uses of the inline data processing.
Cheers,
Gnuplotter

Plot iterations and pseudo-files

2010-01-10T04:59:00.000-08:00

As I promised some time ago, I will discuss some of the new features in gnuplot 4.4. The first one that I would like to show is the concept of iteration in the plot command, and the concept of certain pseudo-files. If you have ever had to script the creation of your plots, you will appreciate these features.

First, let us see what the for loop looks like in the plot command. There are forms of it: once one can loop through integers, while in the other case, one can step through a string of words. So, the iteration either looks like

plot for [var = start : end {:increment}]

plot for [var in "some string of words"]

After this introduction, let us see how these can be used in real life. The first example that I will show is that of waterfall plots, i.e., when a couple of curves are plotted on the same graph, and they are shifted vertically. This is common practice, when one has several spectra, and wants to show the effect of some parameter on the spectra.

reset
f(x,a) = exp(-(x-a)*(x-a)/(1+a*0.5))+0.05*rand(0)
title(n) = sprintf("column %d", n)
set table 'iter.dat'
plot [0:20] '+' using (f($1,1)):(f($1,2)):(f($1,3)):(f($1,4)):(f($1,5)):(f($1,6)) w xyerror
unset table

set yrange [0:15]
plot for [i=1:6] 'iter.dat' u 0:(column(i)+2*i) w l lw 1.5 t title(i)

I would like to walk through the script line by line, for there is something unusual in almost each line. So, after re-setting the gnuplot session, we define a function, which will be a Gaussian, whose centre and width is determined by the parameter 'a'. We then plot this function to a file, 'iter.dat', and do it 6 times, and each time with a different parameter, so that the Gaussian is shifted, and becomes broadened. Note, however, that we do this by plotting a special file, '+'. This was introduced in gnuplot 4.4, and the purpose of this special file is that by invoking this, one can use the standard plot modifiers even with functions. I.e., we can specify 'using' for a function. The importance of this is that many plot styles require several columns, and we could not use those plot styles with functions without the '+' pseudo-file. Consider the following example

reset
unset colorbox
unset key
set xrange [0:10]
set cbrange [0:1]
plot '+' using ($1):(sin($1)):(0.5*(1.0+sin($1))) with lines lw 3 lc palette, \
'+' using ($1):(sin($1)+2):($1/10.0) with lines lw 3 lc palette

which produces the following graph

We can thus colour our curve by specifying the colour in the third column of the pseudo-file. Of course, this is only one possibility, and there are many more. If one wants to plot 3D graphs, then the pseudo-file becomes '++', but the concept is the same: the two variables are denoted by $1 and $2, and the function is calculated on the grid determined by the number of samples and the corresponding data range.

Now, back to the iteration loop! We produce 6 columns of data by plotting '+' by invoking a plot style that requires 6 columns. In this case, it is the xyerrorbars. Having created some data, we plot each column, but we call plot only once: the iteration loop does the rest. In each plot, the curve is shifted upwards, and the title is taken from the column number. For specifying the title, we use the function that we defined earlier: it takes an integer, and returns a string. At the end of the day, we have this graph

This was an example, when we plot various columns from the same file. We can also use the iteration to plot different files. When doing so, there are two options available. One is that we simply specify the file names in a string, as below

reset
filenames = "first second third fourth fifth"
plot for [file in filenames] file using 1:2 with lines

which will plot files 'first', 'second', 'third', 'fourth', and 'fifth'. At this point, note that 'file' is a string, i.e., we can manipulate it as a string. E.g., if we wanted to, instead of 'first', 'second', etc., plot 'first.dat', 'second.dat', and so on, we would do this as

reset
filenames = "first second third fourth fifth"
plot for [file in filenames] file."dat" using 1:2 with lines

The second option is, if the data files are numbered, e.g., if we have 'file_1.dat', 'file_2.dat', and so on, we can use the iteration over integers as follows

reset
filename(n) = sprintf("file_%d", n)
plot for [i=1:10] filename(i) using 1:2 with lines

which will plot the second column versus the first column of 'file_1.dat' through 'file_10.dat'.

Defining some new plot styles

2009-12-03T12:15:00.000-08:00

I have long wanted to write a blog post on this subject, but somehow, I never got the time to do it. But the time has come, and I will do it now. One of the plot styles that I actually like in origin (no matter what I think of the software in overall) is the one where the circumference of a circle is drawn in a colour different to that of the body of the circle, as in this graph

The only glitch is that this is not supported in gnuplot. Well, this is not completely true, because one can use a new prologue for the postscript output, but that works for postscript only, nothing else. What should we do then? There are two options: one has already been discussed, in connection with the bubble charts, sometime in early September. Nevertheless, that method is feasible only when the points do not overlap. What we did there was to plot the data twice (or however many times). But the problem is that those are actually two different plots, so if the adjacent circles overlap, we will not have what we want. Just try it, if you are still not convinced! The second option is that we keep reading. So, the question is then, how could we trick gnuplot into thinking that we have only one plot, not two. Well, the short answer is that we plot only one plot, not two, and we then make sure that points in the plot are coloured properly. I think it is time to show my script, that should make everything clear.

reset
set table 'two.dat'
plot [0:10] sin(x)+0.1*rand(0), cos(x)+0.1*rand(0)
unset table

f(x) = cos(x*pi)
set cbrange [-1:2]; unset colorbox; set border back

set palette defined (-1 "#000000", 1 "#ff0000", 2 "#0000ff")

plot "< gawk '{print $0; print $0}' two.dat" using 1:2:(2+f($0)/5.0):(f($0+1)) index 0 \
with points pt 7 ps var lt palette title 'red', \
'' using 1:2:(2+f($0)/5.0):(f($0+1)*2) index 1 \
with points pt 7 ps var lt palette title 'blue'

As always, we generate some data. Note, that we plot sin(x) and cos(x), which will be in the same data file, but in two separate data blocks. This will become important later. What we need to know about data blocks is that in the data file they are separated by pairs of blank lines, and that we can ask gnuplot to plot only certain data blocks. We indicate our choice with the use of the index keyword. If you want to know more about this, issue the

?index

command, and look at the data file 'two.dat'!

OK, so we have some data, if a bit obscure at the moment. Next, we define a funny function. If you look closely, you will realise that f(x) takes the value of -1 for odd, and 1 for even numbers. In principle, any function should do here, but one has got to be a bit careful: those functions that take 1 or -1 at isolated points might not work. If you do not want to dive into the details of this, take my word for it that f(x) defined above will just be perfect for our purposes.

The next step is the definition of a colour palette with the colour range. (We take off the colour box, too, and set the border to back, but these are irrelevant nuances.) We use the palette for colouring our graph: in the first curve, every other point will be black and red, while in the second curve, every other point will be black and blue. Thus, I have already divulged my trick: we have to duplicate our data set in a way that every line is copied at its position, so, in 'two.dat' we will have something like this

...
 9.29293 -0.93128  i
 9.29293 -0.93128  i
 9.39394 -0.978266  i
 9.39394 -0.978266  i
 9.49495 -0.923256  i
 9.49495 -0.923256  i
...

For this we use an external gawk script, but it is really just one line, it could not be any simpler. Once we duplicated our data, we plot it, and colour every second point as black, and every second point as red. Also note that we use the index keyword to choose the data block that we need. If you have only one data block, you can skip this. But not only do we colour the points differently, but we also resize them: after all, if all had the same size, we would not see the ones that are plotted first. For all this machinery, we use the fact that when a 2D plot is given as 4 columns, the size of the points can be determined by the third column, while the fourth column can be assigned to take care of the colour. In this case, the colour is taken from the palette that we defined above. The general form of such a plot is

plot 'foo' using 1:2:3:4 with points pt 7 ps var lt palette

where 'var' signifies the variable point size (ps), while palette is the colour. Note that we use our snappy f(x) function to choose the size and the colour of every second point, simply by counting the ordinal number of the particular data record, $0: For even numbers, we set the size 2+1/5, and the colour to black (colour value of -1), while for odd numbers, the size is 2-1/5, and the colour is red (colour value of 0). If you have a single data, the whole script would be simply

f(x) = cos(x*pi)
set cbrange [-1:1]; unset colorbox; set border back

set palette defined (-1 "#000000", 1 "#ff0000")

plot "< gawk '{print $0; print $0}' two.dat" using 1:2:(2+f($0)/5.0):(f($0+1)) \
with points pt 7 ps var lt palette title 'red'

Well, this is it for today. Next time I will try to show some of the new stuff in gnuplot 4.4. Cheers,
Zoltán

New version of gnuplot is released!

2009-12-01T12:17:00.000-08:00

I have been waiting for this for a long time, but at long last, it has happened! The new version of gnuplot has been released with the designation 4.4. You can download the binary or the source code from the sourceforge repository. In the future, I will discuss the new things, and show what can be done with them.
Cheers,
Zoltán

Restricting fit parameters

2009-12-01T11:07:00.000-08:00

Chris asked an interesting question today, namely, how one can restrict the fit range in gnuplot. What he meant by that was not the range of the data points (that is really easy, the syntax is the same as for plot), but the range of fit parameters. In some cases, it is a quite reasonable requirement, because we might know from somewhere that certain parameter values just do not make any sense. As it turns out, it is rather easy to achieve this in gnuplot. All we have to do is to come up with a function that restricts its values in the desired range.

After this interlude, let us see an example! We will create some data with the following gnuplot script:

reset
a=1.0; b=1.0; c=1.0
f(x) = a*exp(-(x-b)*(x-b)/c/c)
set table 'restrict.dat'
plot [-2:4] f(x)+0.1*(rand(0)-0.5)
unset table

We take a Gaussian, with some noise added to it. Naturally, we would like to fit a Gaussian to this data, and in particular, f(x). But what, if our model is such that 'a' must be in the range [1.1:2], 'b' must be in the range [0.1:0.9], and 'c' must be in the range [0.5:1.5]? We just use in our fit, instead of f(x), another function, g(x), say, of the form

g(x) = A(a)*exp(-(x-B(b))*(x-B(b))/C(c)/C(c))

where A(a), B(b), and C(c) take care of our restrictions. These functions are somewhat arbitrary, but for better or worse, I will take the following three arcus tangents

# Restrict a to the range of [1.1:2]
A(x) = (2-1.1)/pi*(atan(x)+pi/2)+1.1

# Restrict b to the range of [0.1:0.9]
B(x) = (0.9-0.1)/pi*(atan(x)+pi/2)+0.1

# Restrict c to the range of [0.5:1.5]
C(x) = (1.5-0.5)/pi*(atan(x)+pi/2)+0.5

which would look like this

The point here is that as x runs from negative infinity to positive infinity, A(x) runs between 1.1, and 2.0, and likewise for B(x), and C(x). Then the fit goes on as it would normally. Our script is, thus, the following in its full glory:

# Restrict a to the range of [1.1:2]
A(x) = (2-1.1)/pi*(atan(x)+pi/2)+1.1

# Restrict b to the range of [0.1:0.9]
B(x) = (0.9-0.1)/pi*(atan(x)+pi/2)+0.1

# Restrict c to the range of [0.5:1.5]
C(x) = (1.5-0.5)/pi*(atan(x)+pi/2)+0.5

a=0.0
b=0.5
c=0.9
fit f(x) 'restrict.dat' via a, b, c

g(x) = A(aa)*exp(-(x-B(bb))*(x-B(bb))/C(cc)/C(cc))
aa=1.5
bb=0.5
cc=0.9
fit g(x) 'restrict.dat' via aa, bb, cc

plot f(x), g(x), 'restrict.dat' w p pt 6

and it produces the following graph:

At this point, we should not forget, that what we are interested in is not the value of 'aa', 'bb', or 'cc', but the value of 'a', 'b', and 'c'. This means that what we have to take is
A(aa), B(bb), and C(cc). If you print the values of 'a', 'b', and 'c' from the fit to f(x), the value of 'aa', 'bb', and 'cc', and the value of A(aa), B(bb), and C(cc), we get the following results

gnuplot> pr a, b, c
0.984221984191135 0.996600824504231 1.00765240463672

gnuplot> pr aa, bb, cc
-3442408.91578921 1443864.45093385 -0.201236322474146

gnuplot> pr A(aa), B(bb), C(cc)
1.10000008322047 0.899999823634477 0.936788734177637

Obviously, the second print does not make too much sense, we have to compare the last one to the first one. We can here see that we got values in the ranges [1.1:2], [0.1:0.9], and [0.5:1.5], as we wanted to.

Broken histograms

2009-11-29T10:23:00.000-08:00

Sometimes, a histogram is just a bit awkward, for the simple reason that one or two values are extremely high compared to the rest of the graph. In the case of a standard graph, we would use a broken axis to bring all points to the same order of magnitude. We can play the same trick with histograms, in fact, it is, in some sense, even simpler, than the broken axes. All we have to do is to plot a thick line at the proper position in the proper colour. This is the graph that we are going to make today

Our data file, brokenhist.dat, is as follows

"Jan" 1 2
"Feb" 44 4
"Mar" 3 1
"Apr" 2 25
"May" 4 5
"June" 2 1

and here is our script:

reset
blue = "#babaff"
set xrange [-0.5:5.5]
set yrange [0:11]
set isosample 2, 100
set table 'brokenhist_b.dat'
splot 1-exp(-y/2.0)
unset table

unset colorbox
set border 1
set xtics rotate by 45 nomirror offset 0, -1
unset ytics
f(x) = (x < 6 ? x : (x < 30 ? x-17 : x-35) )
g(x) = (x < 6 ? 1/0 : 6)
set boxwidth 0.85
set style fill solid 0.8 border -1
set style data histograms
set palette defined (0 "#ffff   ff", 1 "#babaff")
plot 'brokenhist_b.dat' w ima,\
'brokenhist.dat' u (f($2)) t 'Red bars',\
'' u (f($3)) lc rgb "#00bb00" t 'Green bars', \
'' u 0:(f($2)):2 w labels center offset 0,0.5 t '',\
'' u ($0+0.25):(f($3)):3 w labels center offset 0,0.5 t '',\
'' u 0:(-1):xticlabel(1) w l t '', \
'' u ($0+0.12):(g($3)+0.12):(0.25):(0.25) w vectors lt -1 lc rgb blue lw 5 nohead t '', \
'' u ($0-0.12):(g($2)+0.12):(0.25):(0.25) w vectors lt -1 lc rgb blue lw 5 nohead t ''

OK, so let us look at the code! The first couple of lines are required only, if you want to have some posh background. Likewise, you can drop the 'unset colorbox' line, when you have a white background. We set only the bottom axis, which means that we have to unset the ytics and set to xtics to nomirror. Then we have two helper functions. The definitions of these depend on where you want to have the break point in the histogram. In this particular case, I took 6, but it is arbitrary.

In the next step, we set the properties of the histogram, like the width of the columns, the fill style, and the data style. We also define a palette, but this is needed for the background only. For white background, you can skip this step. You can also skip the first plot, because that is nothing but our fancy background.

The actual plotting begins after this. We plot the two sets of columns, and also plot the data file with labels. The labels are placed at the top of each column (this is why we could do away with the yaxis.) We also 'plot' the axis labels, and finally, plot the two break points. Note that the plotting of the break points is automatic, once we have the definitions of the two helper functions. If you want to have a steeper cut, you could

'' u ($0+0.12):(g($3)+0.12):(0.25):(0.5) w vectors lt -1 lc rgb blue lw 5 nohead t ''

e.g., which stretches the vectors in the vertical direction. Otherwise, we have finished the plot, there is nothing else to do.
I should point out here that, in order to have a seamless cut, we have to use a colour for the vectors that is identical to the background at that particular point. This implies that we could not have a gradient at y=6. The background colour is virtually constant at y=6 (c.f. the definition of 'brokenhist_b.dat'. While it would not be impossible to implement a cut over a gradient, I believe, it is probably not worth the trouble it involves.

Update

2009-11-22T11:04:00.000-08:00

I have had some time, so I moved all recent posts to their permanent place on
my web page. I have "sexed up" the homepage a bit, so, hopefully, browsing will be a tad easier. Let me know, if there are any problems! (I know that there is a small glitch with the cascaded style sheets on IE6. IE8 should work without problems. Firefox 3.5 is also OK.) There is a zipped version of the complete site, if you want to read it off-line.
Comments should still be posted here, please!

Cheers,
Zoltán

Some basic statistics with gnuplot

2009-11-22T00:54:00.000-08:00

In my previous post, I mentioned a patch that you can compile into gnuplot, and that should make plots with some statistical properties a bit easier. Now, the problem with that patch is that, if you don't want to, or can't take the trouble of compiling gnuplot for yourself, it is no use. However, for most things contained in the patch, there is a work-around that should function properly in gnuplot 4.2. I will discuss those now.

The first thing I did with the statistical patch was to plot the mean, minimum and maximum of a data set. This can easily be done in the following way.

reset
# Produce some dummy data
set sample 200
set table 'stats2.dat'
plot [0:10] 0.5+rand(0)
unset table

set yrange [0:2]
unset key

# Retrieve statistical properties
plot 'stats2.dat' u 1:2
min_y = GPVAL_DATA_Y_MIN
max_y = GPVAL_DATA_Y_MAX

f(x) = mean_y
fit f(x) 'stats2.dat' u 1:2 via mean_y

# Plotting the minimum and maximum ranges with a shaded background
set label 1 gprintf("Minimum = %g", min_y) at 2, min_y-0.2
set label 2 gprintf("Maximum = %g", max_y) at 2, max_y+0.2
set label 3 gprintf("Mean = %g", mean_y) at 2, max_y+0.35
plot min_y with filledcurves y1=mean_y lt 1 lc rgb "#bbbbdd", \
max_y with filledcurves y1=mean_y lt 1 lc rgb "#bbddbb", \
'stats2.dat' u 1:2 w p pt 7 lt 1 ps 1

At the beginning of our script, we just produce some dummy data, and call a dummy plot. This plot does nothing but fills in the values of the minimum and maximum of the data set. Then we fit a constant function. You can convince yourself that this returns the average of the data set.

In the plotting section, we produce three labels that tell us something about the data set, and plot the data range with shaded region. Easy enough, and in just a couple of lines, we created this figure

Now, what should we do, if we were to calculate the standard deviation. Well, we know how to calculate the average, so we will use that. Here is the script:

reset
set sample 200
set table 'stats2.dat'
plot [0:10] 0.5+rand(0)
unset table

set yrange [0:2]
unset key
f(x) = mean_y
fit f(x) 'stats2.dat' u 1:2 via mean_y

stddev_y = sqrt(FIT_WSSR / (FIT_NDF + 1 ))

# Plotting the range of standard deviation with a shaded background
set label 1 gprintf("Mean = %g", mean_y) at 2, min_y-0.2
set label 2 gprintf("Standard deviation = %g", stddev_y) at 2, min_y-0.35
plot mean_y-stddev_y with filledcurves y1=mean_y lt 1 lc rgb "#bbbbdd", \
mean_y+stddev_y with filledcurves y1=mean_y lt 1 lc rgb "#bbbbdd", \
mean_y w l lt 3, 'stats2.dat' u 1:2 w p pt 7 lt 1 ps 1

What we utilise here is the fact that the fit function also sets a couple of variables. One of them is the sum of the residuals, which is called FIT_WSSR, while another is the number of degrees of freedom, FIT_NDF. However, we know that the number of degrees of freedoms is one less, than the number of data points, for we fit a function with a single parameter. Therefore, if we take the square root of the sum of residuals divided by the number of degrees of freedom plus one, we get the standard deviation. The rest of the plot is trivial, and this script results in the following graph:

Incidentally, this can also be used for removing points that are very far from the mean. The following script takes out those data that are more than one standard deviation away from the mean.

reset
set sample 200
set table 'stats2.dat'
plot [0:10] 0.5+rand(0)
unset table

set yrange [0:2]
unset key
f(x) = mean_y
fit f(x) 'stats2.dat' u 1:2 via mean_y

stddev_y = sqrt(FIT_WSSR / (FIT_NDF + 1 ))

# Removing points based on the standard deviation
set label 1 gprintf("Mean = %g", mean_y) at 2, min_y-0.15
set label 2 gprintf("Sigma = %g", stddev_y) at 2, min_y-0.3
plot mean_y w l lt 3, mean_y+stddev_y w l lt 3, mean_y-stddev_y w l lt 3, \
'stats2.dat' u 1:(abs($2-mean_y) < stddev_y ? $2 : 1/0) w p pt 7 lt 1 ps 1

with the corresponding figure

Only the last line is relevant: we use the ternary operator to decide whether we want to keep the point: if the deviation from the mean is less, than the standard deviation, we hold on to our data, otherwise, we replace it by 1/0, which is undefined, and gnuplot quietly ignores it. If you want to learn more about the working of the ternary operator, check out my post on the plotting of an inequality.

We have, thus, already found a solution for two of the problems addressed in the patch. What about the third one, adding arrows to the plot at the position of the minimum or maximum, say. We can do that, too. Here is the script:

reset
set sample 50
set table 'stats1.dat'
plot [0:10] 0.5+rand(0)
unset table

set yrange [0:2]
unset key
plot 'stats1.dat' u 1:2
min_y = GPVAL_DATA_Y_MIN
max_y = GPVAL_DATA_Y_MAX

plot 'stats1.dat' u ($2 == min_y ? $2 : 1/0):1
min_pos_x = GPVAL_DATA_Y_MIN
plot 'stats1.dat' u ($2 == max_y ? $2 : 1/0):1
max_pos_x = GPVAL_DATA_Y_MAX

# Automatically adding an arrow at a position that depends on the min/max
set arrow 1 from min_pos_x, min_y-0.2 to min_pos_x, min_y-0.02 lw 0.5
set arrow 2 from max_pos_x, max_y+0.2 to max_pos_x, max_y+0.02 lw 0.5
set label 1 'Minimum' at min_pos_x, min_y-0.3 centre
set label 2 'Maximum' at max_pos_x, max_y+0.3 centre
plot 'stats1.dat' u 1:2 w p pt 6

First, we retrieve the values of the minimum and the maximum by using a dummy plot. Having done that, we retrieve the positions of the minimum and maximum, by calling a dummy plot on the columns

plot 'stats1.dat' u ($2 == min_y ? $2 : 1/0):1

What this line does is substitute min_y, when the second column (whose minimum we extracted before) is equal to the minimum, and an undefined value, 1/0, otherwise. The minimum of this plot is nothing, but the x position of the first minimum. Likewise, had we assigned

min_pos_x = GPVAL_DATA_Y_MAX

that would have given the position of the last minimum of the data file. Obviously, these distinctions make sense only, if there are more than one minimum or maximum. Knowing the x and y positions of the minimum and maximum, we can easily set the arrows. We, thus, have the following figure

Adding labels showing the value should not be a problem now.

Well, this is for today. Till next time!

Patching gnuplot

2009-11-09T12:21:00.000-08:00

One of the major advantages of open-source code is that if you would like to add some new features, you can easily do that. This applies to gnuplot, too, in fact, doing that does not require anything special. I have been quite inactive on this blog recently, and the reason is that Philipp Janert and I have been working on a patch to gnuplot.

The steps of patching gnuplot are described on gnuplot's main web page. There are a number of patches uploaded to gnuplot's patch tracker, on which quite a few new features, still in the development phase, are published. It is really worthwhile to try them out, first, to provide feedback as to what is useful and what is not, and second, to help the developers to find bugs and other glitches, like what the syntax of a command should be and so on.

Our patch is related to an old debate as to what gnuplot really is. At many a place, you will find the statement that "gnuplot is a plotting utility, not a statistical analysis package". I have nothing against this statement, however, when saying so, we have to tell what we mean by plotting. So, is plotting just placing a thousand dots at positions that represent our data? Or do we want more? E.g., throwing out data points that are unreasonably far from the mean. Or showing the mean, and the standard deviation? Or calling the reader's attentional to some special points, like the minimum or the maximum in a data set? And many similar things. I believe, plotting requires much more, than just showing the measurement data: a plot makes sense only, if we can point out what is to be pointed out. By the way, fitting falls into this category, and fitting has been an integral part of gnuplot for ages. The point being that the original statement (gnuplot is a plotting utility, not a statistical analysis package) has been wrong for a long time.

The patch that I mentioned above was announced yesterday on the gnuplot development mailing list and you can find the patch for the source and the
documentation on patch tracker. I have put a couple of examples on my gnuplot web site under patch. You can also find the full documentation.

I would like to ask you, if you feel crafty and you can, download the patch, and try it, and let us know whether you find it useful, what else, do you think, we could do with it and so on. It would really help the development. Once the patch makes it to the main code, I will discuss various option on these pages.

Just to wet your appetite, here is a figure that you could very easily make with the new patch. (You can find the code on my web site.)

Many cheers,
Zoltán

A long shadow is cast on the plot...

2009-10-16T12:34:00.000-07:00

A couple of days ago, Cedric asked the question how we could add a shadow to a pm3d surface plot. In some sense, the problem turned out to be easier, than I had expected, but I am not sure that I did what he actually meant. Anyway, we could use it as our working hypothesis, and refine it, if necessary.

The idea of putting "phong" on the surface was discussed ages ago, and I won't re-open that question here. For the shadow, we will just replot our surface (defined as z(x,y)) in a little bit strange way: instead of letting the x and y run through their corresponding ranges, we will restrict y to be equal to the maximum of the yrange. By doing so, we get this

Here is our (short) script:

reset
unset colorbox; unset key
set iso 2, 50
set parametric; set urange [0:1]; set vrange [0:1.2]
set table 'tb.tab'
splot u, v, 1
unset table
unset parametric
set iso 100, 100
set xrange [-3:3]
set yrange [-3:3]
set zrange [0:1.2]
set table 't.tab'
splot exp(-x*x-y*y)+0.1*rand(0)
unset table
set ticslevel 0
set pm3d
set cbrange [0:4]
f(x,y,z,a,b,s) = z*(exp(-(x-a)*(x-a)/s-(y-b)*(y-b)/s)/3.0+0.66)
set palette defined (0 "#ff2222", 1 "#ffeeee", 2 "#aaaaaa", 3 "#2222ff", 4 "#8888ff")
splot 'tb.tab' u ($1*6.0-3.0):(3):2:($2+3.0) w pm3d, \
'' u (-3):($1*6.0-3.0):2:($2+3.0) w pm3d, \
'' u ($1*6.0-3.0):($2*6.0-3.0):(0):(2.1) w pm3d, \
't.tab' u ($1+1.0):(3):($3*0.7):(2.1) w pm3d,\
'' u 1:2:3:(f($1,$2,$3,-0.5,-0.5,.8)) w pm3d

First, we create the data that will be our background, then some dummy data, which in this particular case will be a noise Gaussian function in 2D. Then we move the surface to the bottom of the zrange, by setting the ticslevel to 0. The next step is the definition of the cbrange. We need to "overdefine" this, i.e., our cbrange is much larger, than the actual data range. The reason for this is that in this way, we can use the same colour palette, and we do not have resort to multiplot. The basic problem with multiplot is that we since we want to plot onto the same graph, we would have to re-set the border, tics, labels, and so on. By using the same plot, and same palette, we can avoid all this hassle. The price we pay is that our palette will be a bit more complicated, but we can live with that. Now, we have to define our palette, which will change between red and almost white for [0:1], and between blue, and almost white for [3:4]. We also defined a value at 2, which we will use for colouring the shadow, but the two main ranges are [0:1], and [3:4]. Note that we define disjoint ranges, thereby not walking into a trouble with the end points.
Finally, we plot our background, and function. Pay attention to plots like

splot 'tb.tab' u ($1*6.0-3.0):(3):2:($2+3.0) w pm3d

This will plot a plane between x [-3:3], and z [0:1] at y=3, with a colour given by the value of ($2+3.0). This is nothing but pushing all z values in the plot into the [3:4] colour range.

At the very end, we plot our function's shadow, namely,

't.tab' u ($1+1.0):(3):($3*0.7):(2.1) w pm3d

where the x values are shifted by 1.0 (so as to give the impression that the surface is lit from the (-1:-1) direction), the y values are all restricted to 3, which is the maximum of the yrange, and the z values are multiplied by 0.7, again for the same reason. Colouring is done by using one single value, 2.1. The very last step is to plot the surface itself, using the colouring given by f(...). If you want to have another direction for the lighting, or have a tighter focus, you should change the parameters in this function.

The shiny histograms, again

2009-10-13T13:12:00.000-07:00

Do you remember those shiny histograms that we discussed some long, long time ago, at the beginning of the summer? And do you also remember how much of a hassle it was to create them, since we relied on an external gawk script, and we had to build our rectangles, one by one? And have you ever wondered what all that fuss was about and whether there was an easier solution? A one-liner, perhaps? If the answer to these questions lies in the affirmative, go no further! We will discuss a method of making those histograms, without an external script, only with legal gnuplot commands, and in 5 lines. I understand that 5 lines is just 4 lines longer, than one would expect from a one-liner, but on the other hand, three out of those 5 lines are equivalent, so I don't feel so bad about this any more:)
OK, so here is our figure

here is our data file, which we will call 'hist.dat'

and here are our scripts, 'hist.gnu',

reset
unset key; set xtics nomirror; set ytics nomirror; set border front;
div=1.1; bw = 0.9; h=1.0; BW=0.9; wd=10; LIMIT=255-wd; white = 0
red = "#080000"; green = "#000800"; blue = "#000008"
set auto x
set yrange [0:11]
set style data histogram
set style histogram cluster gap 1
set style fill solid
set boxwidth bw
set multiplot
plot 'hist.dat' u 1 lc rgb red, '' u 2 lc rgb green, '' u 3 lc rgb blue
unset border; set xtics format " "; set ytics format " "; set ylabel " "
call 'hist_r.gnu'
unset multiplot

and 'hist_r.gnu'

bw=BW*cos(white/LIMIT*pi/2.0); set boxwidth bw; white=white+wd
red = sprintf("#%02X%02X%02X", 128+white/2, white, white);
green = sprintf("#%02X%02X%02X", white, 128+white/2, white);
blue = sprintf("#%02X%02X%02X", white, white, 128+white/2);
rep
if(white<LIMIT) reread

Then let us see what is happening here! At the beginning, we define various variables, most notably, white, red, green, and blue. The rest up to the first plot command is nothing but setting up the figure: we define the range, tell gnuplot to treat our data as histogram, set the width of the bars, and finally, set multiplot.
There is nothing exciting in the first plot, except, that we specify the colour of the bars as

plot 'hist.dat' u 1 lc rgb red, '' u 2 lc rgb green, '' u 3 lc rgb blue

The strings red, green, and blue were defined at the beginning of our first script, thus, we learn here that gnuplot will accept any defined (and valid) string as the specifier of the colour. After our first plot, we unset the border, re-set the format of the xtic and ytic to empty, and do likewise with the ylabel. Should there be an xlabel, we should have to do the same there. Having done this, we call our second script, which we will dissect now. This is really nothing but a 'for' loop, that we have discussed a couple of times before. In fact, quite a few times. The first two commands re-set the widths of the bars in the next plot. Note that I set the width in such a way that it would draw the outline of a circle as we step through the values of white. This is what we increment next, mind you!

The next three lines are basically identical: we re-define the colours, using a sprintf command in each step. If you recall how the RGB colours are defined, we have to create a string that looks like

#00FF00

say. This is what our sprintf command will do, returning a string of this form that depends on the value of 'white'. There are some small nuances in the corresponding colour channel of red, green, and blue, respectively, namely, that the base colour for red was

#080000

and we want to linearly interpolate between this colour, and white,

#FFFFFF

so, we have to apply the relevant linear function, but there is nothing beyond this. Obviously, if you are unhappy with the colour scheme that I have (I know full well that these are not the best colours...), this is the place where you would have to tamper with the script. When we are done with re-defining the widths, and the colours, we simply replot our histogram, and do that, as long as the value of white is smaller, than the limit that we set at the beginning. In this particular case, 245. In most cases, we do not need this many plots, by the way! For a raster plot, 10-12 steps, for a vector format, something like 15-16 steps should be more than enough.

At the end, we shouldn't forget about unsetting the multiplot. The only difficulty that I see with this figure is that it is not so straightforward to use a key. However, it is not terribly hard to come up with a solution for this problem: all we have to do is to put three vertical labels on the top of the first, second, and third column, indicating what they represent.

How useful was this plot?

Putting figures into perspective

2009-10-12T12:13:00.000-07:00

As I promised yesterday, this time, we will try to skew our figure, as if it was in 3D, and we were looking at it from an angle. I must admit that this is something that one would not put in a publication, unless it is a poster or presentation, perhaps. In those cases, however, it might lend a refreshing new, errr.., perspective to our data. Before diving into the script, here is the figure, so you can decide whether you want to read on:)

Then, here is our script that was responsible for the figure

reset
xmin = 0; xmax = 10; zmin = -0.4; zmax = 0.75
set view 60, 30
unset key; unset colorbox
set border 1+16+128+1024
unset ytics; set xtics out nomirror; set ticslevel 0
set yrange [0:-0.1]
set zrange [zmin:zmax]
set grid front; unset grid
set xtics xmin,2,xmax-1
set ztics zmin,0.2,zmax
f(x) = exp(-x/4.0)*sin(x)
c(x) = exp(-(x-xmax/3.0)*(x-xmax/3.0)/1.0)
set xlabel 'Time [a.u.]'
set label 2 'Amplitude [a.u.]' at graph -0.35, 0.3 rotate by 90
set parametric
set iso 3, 3
set urange [xmin:xmax]
set vrange [zmin:zmax]
set table 'perspective1.dat'
splot u, 0, v
unset table
set urange [xmin+0.2:xmax-0.2]
set table 'perspective2.dat'
splot u, 0, f(u)
unset table
unset parametric
set size 1.4, 1.4
set palette defined (0 "#7b6cff", 1 "#eeeeff")
splot 'perspective1.dat' u 1:2:3:(c($1)) w pm3d, \
'perspective2.dat' u ($1+0.05):2:($3-0.02) w l lw 6 lc rgb "#888888", \
'' w l lw 6 lt 1

The trick that we use here is to plot in 3D, and take off all those elements of the figure that we do not actually need. In the beginning, we define a couple of variables, in order to make our life a bit easier. Then we unset the colorbox and the key, and set only those borders that we want to see. If you are interested in how I came up with those numbers in the definition of the border, just issue the command

?border

in the gnuplot prompt. After this, we keep unsetting things, and define our ranges. The next noteworthy command is

set grid front; unset grid

At first, this seems silly, but the reason is real: in the figure we have a background, which would obscure our tic marks. The way out of this problem is to push the tic marks to the front, which is achieved by setting the grid to the front. Since we do not actually need the grid, we unset it, but the setting, its front position, is still there. The tic marks inherit the position of the grid, and thus, they will be in the forefront.

The function definitions are f(x), the function that we want to plot, and c(x), which will determine the colouring of our background. Modify them accordingly.

We set the zlabel by hand, because it might not be possible to turn it by 90 degrees otherwise. (Not all terminals would support it.) Then we plot the background, and the function, both to a file, so that it will be easier to colour them in the next step, where we plot the real thing. Note that in the plot of the background, we use four columns, while there are only 3 in the data file. The fourth one determines the colour as the position on the graph. The colour is given by the function c(x), and the palette, which we defined one line earlier. You should change these two things, if you are not satisfied with the background you get. Finally, we plot the function twice. Once in gray, a bit shifted to the right and down, and for the second time, in red. In this way, we add a shadow to our curve. If you want to improve the shadow, you should look at my post from the 30th August

I should also add that we discussed a vertical skew. In case you want to skew the figure horizontally, all you have got to do is to plot on the x-y plain instead of x-z.

How useful was this plot?

Adding a variable grid to a plot

2009-10-11T05:50:00.000-07:00

Useful as it is, sometimes the standard grid is not the best way to guide the eye when reading off values in a graph. Today I would like to discuss a couple of other options, none of them standard in gnuplot. Fortunately, we haven't got to leave the realm of gnuplot, everything can be done without any external scripts. To make things even more appealing, the solution requires only a couple of lines of code.

Our first attempt will be to "project" the values of a curve to the two axes, as in this figure

The data file, helpergrid.dat, contains only a couple of lines,

and the script reads as follows

reset
unset key
set xrange [0:40]
set yrange [0:7]
set xlabel 'Effort [a.u.]'
set ylabel 'Output [a.u.]'
p 'helpergrid.dat' u 1:2:(0):1 w xerror ps 0 lt 0 lw 0.5 lc rgb "#ff0000", \
'' u 1:2 w i lt 0 lw 0.5 lc rgb "#ff0000", \
'' u 1:2 w lp pt 7 ps 1.5

Drawing the vertical lines is easy, for there is a plotting style in gnuplot, impulse, which does just that. The horizontal ones are not that trivial, but take only one line of code: the idea is that we trick gnuplot into drawing the horizontal lines by something else. In this case, this something else is the error bar, in particular, the xerror bar. This is what happens in the first plotting line, after we set up the figure. The first plot calls our data file, and draws the errors that span from x=0 to x=x_data. We also specify the line type, 0, which will be dashed, the line width, and the line colour. The second plot draws the vertical lines, plotting our data with impulses. Again, we specify the line type, width, and colour, so that it conforms with the previous one. This is necessary, because, when left alone, gnuplot assigns a new line type to each new plot. Finally, we plot the data. This was simple.

I should add here that we can draw a "full" grid based on the data values. If we replace the plotting command by

p 'helpergrid.dat' u (40):2:(0):(40) w xerror ps 0 lt 0 lw 0.5 lc rgb "#ff0000", \
'' u 1:(7) w i lt 0 lw 0.5 lc rgb "#ff0000", \
'' u 1:2 w lp pt 7 ps 1.5

then the horizontal lines will be drawn between 0 and 40, while the vertical ones between 0 and 7. Now, it might happen that you don't know the xrange and yrange, in which case you can use the variables GPVAL_X_MIN, GPVAL_X_MAX, GPVAL_Y_MIN, GPVAL_Y_MAX to specify the plot. Of course, we have to call a dummy plot beforehand. We have discussed this a couple of times. If in doubt, look up the post on how to plot a martix with bars.

Well, we have cleared the first barrier, but what if we wanted to place the numbers where the data points are? This is quite simple: we only have to plot the file twice more, this time with xticlabel and yticlabel. So, here is our script

reset
unset key
set xrange [0:40]; set yrange [0:7]
set xlabel 'Effort [a.u.]'
set ylabel 'Output [a.u.]'
p 'helpergrid.dat' u 1:2:(0):1 w xerror ps 0 lt 0 lw 0.5 lc rgb "#ff0000", \
'' u 1:2 w i lt 0 lw 0.5 lc rgb "#ff0000", \
'' u 1:2 w lp pt 7 ps 1.5, \
'' u 1:(0):xticlabel(1) w p ps 0, '' u (0):2:yticlabel(2) w p ps 0

and here is our figure

Note that in this case, if you want to produce a full grid (wall-to-wall, floor-to-ceiling), we can just set the grid, and forget about the error bars. The script for that case would be

reset
unset key
set xrange [0:40]; set yrange [0:7]
set xlabel 'Effort [a.u.]'
set ylabel 'Output [a.u.]'
set grid lt 0 lw 0.5 lc rgb "#ff0000"
p 'helpergrid.dat' u 1:2 w lp lt 3 pt 7 ps 1.5, \
'' u 1:(0):xticlabel(1) w p ps 0, '' u (0):2:yticlabel(2) w p ps 0

and this would result in this figure.

I think what we covered today was fairly simple, but still useful. Next time I will show how figures can be skewed, yet another way of making the reader perceive a 2D plot in 3D surroundings. Till then!

More on histograms

2009-10-03T13:38:00.000-07:00

Someone asked me whether it would be possible to make a histogram that looks like those that MS Office can create: filled with gradient, casting a shadow, and labelled according to their value. Frankly, I just couldn't admit defeat, and I had to figure out a way. But then, it turned out to be rather simple, so I thought that I would share it here. We are going to make this figure

from this data, called 'msbar.dat'

Max 1.0 1.0 1.2
Min 0.9 0.2 0.95
Avg 0.95 0.5 1.1

I have already discussed a way to beefing up the histograms (It was titled Shiny histograms, or something similar.), but I must say, that method is a bit convoluted. The reason for this is that we used a parametric plot. We can avoid that by plotting to a file first, and then plotting the file as many times as needed. We can't save the trouble of having to read our data file, but this is rather simple. We can do that either by employing a very primitive external script, or writing the script in gnuplot, as we discussed it in, say, the post on the recession graph. Given that we haven't got to process any data, this latter method is probably less desirous. Gnuplot is not for printing lines and the like, after all.

After the introduction, let us see the script!

reset
# First, the gradient for the bars
set xrange [0:1]; set yrange [0:1]; set isosample 2, 200
set table 'msbar_bar.dat'
splot y
unset table

# Then we make the shadow
set isosample 200, 50
set table 'msbar_bar_sh.dat'
splot (1.0-exp((x-1.0)*20.)) #*(1.0-exp((y-1.0)*20.0))
unset table

reset
unset key; unset colorbox
xw = 0.2; set boxwidth xw; sh = 0.03
gf(x) = x*x*x*x #change this, if a tighter gradient is needed
set cbrange [0:7]
set xrange [-0.5:2.5]
set yrange [0:1.4]
set grid ytics lw 0.5 lc rgb "#868686"
set xtics nomirror
set ylabel 'Performance [a.u.]'
set palette defined (0 '#e0e8f5', 1 '#31bd71', 2 '#e0e8f5', 3 '#d99795', 4 '#e0e8f5', 5 '#9ab5e4', 6 "#ffffff", 7 "#a2a2a2")
plot 'msbar.dat' u 0:(0):xticlabel(1) w l, \
'' u ($0-xw):($2+0.1):(stringcolumn(2)) w labels, \
'' u ($0):($3+0.1):(stringcolumn(3)) w labels, \
'' u ($0+xw):($4+0.1):(stringcolumn(4)) w labels, \
'msbar_bar_sh.dat' u ($1*xw-1.5*xw+sh):($2*1.0-sh):($3+6.0) w ima, \
'' u ($1*xw-.5*xw+sh):($2*1.0-sh):($3+6.0) w ima, \
'' u ($1*xw+.5*xw+sh):($2*1.2-sh):($3+6.0) w ima, \
'' u ($1*xw-1.5*xw+sh+1):($2*0.9-sh):($3+6.0) w ima, \
'' u ($1*xw-.5*xw+sh+1):($2*0.2-sh):($3+6.0) w ima, \
'' u ($1*xw+.5*xw+sh+1):($2*0.95-sh):($3+6.0) w ima, \
'' u ($1*xw-1.5*xw+sh+2):($2*0.95-sh):($3+6.0) w ima, \
'' u ($1*xw-.5*xw+sh+2):($2*0.5-sh):($3+6.0) w ima, \
'' u ($1*xw+.5*xw+sh+2):($2*1.1-sh):($3+6.0) w ima, \
'msbar_bar.dat' u ($1*xw-1.5*xw):($2*1.0):(gf($3)) w ima, \
'' u ($1*xw-.5*xw):($2*1.0):(gf($3)+2.0) w ima, \
'' u ($1*xw+.5*xw):($2*1.2):(gf($3)+4.0) w ima, \
'' u ($1*xw-1.5*xw+1):($2*0.9):(gf($3)) w ima, \
'' u ($1*xw-.5*xw+1):($2*0.2):(gf($3)+2.0) w ima, \
'' u ($1*xw+.5*xw+1):($2*0.95):(gf($3)+4.0) w ima, \
'' u ($1*xw-1.5*xw+2):($2*0.95):(gf($3)) w ima, \
'' u ($1*xw-.5*xw+2):($2*0.5):(gf($3)+2.0) w ima, \
'' u ($1*xw+.5*xw+2):($2*1.1):(gf($3)+4.0) w ima, \
'msbar.dat' u ($0-xw):2 w boxes lt -1 lw 0.5 lc rgb "#4f81bd", \
'' u ($0):3 w boxes lt -1 lw 0.5 lc rgb "#4f81bd", \
'' u ($0+xw):4 w boxes lt -1 lw 0.5 lc rgb "#4f81bd"

If you look at the graph, we have three different gradients, and the shadow. That makes four colour schemes altogether. Since we would like to save the difficulties associated with multiplot, we have to cram all those colour schemes into one colour range. More on this a bit later.

So, first we plot the gradient that will fill the bars, and then the "surface" that will represent our shadow. Then we set up our figure: we take off the keys, and unset the colour box. We also specify the width of our bars 'xw', and set the box width accordingly. This latter step is needed, because we want to have a border to the bars. The shadow shift, 'sh' is also defined here. In the next line, we set our colour range, in this particular case, between [0:7]. As we have pointed out, we need 4 colour ranges, and they should not overlap, therefore, we need 3 gaps between these ranges. For the sake of simplicity, we define the gap to be 1, and all the ranges to be one. This is why we end up with a total colour range of [0:7]. If you have more, or less bars to plot, you can re-define this range accordingly. The next couple of steps are trivial, up to the definition of the colour schemes. We want to have only simply gradients, thus, it is enough to define the colours at the end points. If you are unhappy with the colouring of your bars, you should change the colours here.

Having set up the figure, we can plot the data. First, we plot the labels for the xtics, and the values of the bars. Next comes the shadow. It might be a bit of an overkill for this figure, so you can skip all lines up to 'msbar_bar.dat'. Note that we simply go through all the data points in our file, and shift the shadow in each step. The height of the shadow is given by the product of the second column (which takes values between 0 and 1) of 'msbar_bar_sh.dat' and the particular data point in 'msbar.dat'. We also add a small downwards and rightwards shift to the shadows, lest they should be covered by the bars. The most important point, however, is that we plot the shadow as

'' u ($1*xw-.5*xw+sh):($2*1.0-sh):($3+6.0) w ima

i.e., the third column is shifted by 6. We do this, so that the shadow is pushed into the [6:7] range, where the shadow's colour scheme is defined. Plotting the bars happens in a similar fashion, the only difference is that we do not add 'sh' to the columns, and we push the values into the range of the appropriate colour range. At the very end, we plot the data file once more, this time with boxes, so that the bars can have a border to them.

If you want to use this plot many times, it might be worthwhile to implement it in a script in the language of your choice. The only thing required is printing lines and numbers. In pseudocode, it would look something like this:

for i
  for j
     d = read(datafile(i,j))
     print "'msbar_bar_sh.dat' u ($1*xw-(%d-1.5)*xw+sh):($2*%d-sh):($3+6.0) w ima, \", i, d
     print "'msbar_bar.dat' u ($1*xw-(%d-1.5)*xw):($2*%d):($3+%d) w ima, \", i, d, i
  end
end

Now, suppose that you want to add a legend to the figure. The usual way, setting the key, will not work here, for obvious reasons. However, we can easily add that to the figure. We just have got to find some space for it. Since there is no space left on the figure, we have to make some: we will use multiplot, and specify the size of the main figure to be less, than 1. Therefore, adding the following lines to our gnu script should produce the legend.


set multiplot
set size 0.8, 1
... (The main plot, identical to our original plot)
set origin 0.75, 0
unset border; unset xtics; unset ytics; unset ylabel; unset xlabel
set label 1 at first -0.15, 1.35 "Flatland"
set label 2 at first -0.15, 1.15 "Curvedland"
set label 3 at first -0.15, 0.95 "Land"
plot 'msbar_bar.dat' u ($1*xw-0.4):($2/10.0+1.3):(gf($3)) w ima, \
'msbar_bar.dat' u ($1*xw-0.4):($2/10.0+1.1):(gf($3)+2) w ima, \
'msbar_bar.dat' u ($1*xw-0.4):($2/10.0+0.9):(gf($3)+4) w ima

unset multiplot

The script above results in the figure below:

The turning of the histogram

2009-10-02T11:45:00.000-07:00

I have been off for some time now, and I thought that it would be high time to post something on gnuplot. I have chosen an easy subject, but one that could turn out to be useful. At least, I have seen people searching for a solution to this problem with the histograms. If you needed it in the past, you have probably realised that gnuplot can make only vertical histograms. But sometimes, this is not the best way to represent data, and a horizontal version would be much better. For one thing, the horizontal one might take up much less space. If you want to learn how to produce the figure below, keep reading!

Making the graphs will require some work, but everything is quite straightforward. What we should note, however, is that we will still produce an upright image that we have to rotate later. If you set the terminal to postscript, you will probably use the figure in LaTeX, where you simply have to tell the compiler to rotate the figure by 90 degrees. If you set the terminal to a raster format, or print the screen, you can rotate the image in many applications, but if you want to adhere to the command line, you can, at least, under Linux, use the convert command as

convert -rotate 90 figure_in.png figure_out.png

The figure above was produced based on the following data file

(We have seen these data many times, you should know it by heart by now!) After this interlude, let us see the code! I will discuss it afterwards.

reset
set key at graph 0.24, 0.85 horizontal samplen 0.1
set style data histogram
set style histogram cluster gap 1
set style fill solid border -1
set boxwidth 0.8
set xtic rotate by 90 scale 0
unset ytics
set y2tics rotate by 90
set yrange [0:0.35]; set xrange [-0.5:6.5]
set y2label 'Output' offset -2.5
set xlabel ' '
set size 0.6, 1
set label 1 'Year' at graph 0.5, -0.1 centre rotate by 180
set label 2 'Nowhere' at graph 0.09, 0.85 left rotate by 90
set label 3 'Everywhere' at graph 0.2, 0.85 left rotate by 90
p 'pie.dat' u 2 title ' ', '' u ($2/2.0+rand(0)/10.0) title ' ', '' u 0:(0):xticlabel(1) w l title ''

The first line after the reset is required, because we have to make our key ourselves. We simply specify the coordinates and that that we want to have a horizontal key (i.e., one in which the keys are placed to the right of the previous one), with a length of 0.1. If you want smaller or larger space between the keys, you can modify this, by adding the spacing flag to the set command. You can see the exact syntax by issuing ?key in the gnuplot prompt.
In the next four lines, we set up our histogram. The content of these lines might depend on what exactly we intend to plot. For more on this, check out the gnuplot demo page! We then rotate the xtics by 90 degrees. We do this, because the whole image will be rotated, and by rotating the xtics, we make sure that those will be horizontal at the end. For the very same reason, we also unset the ytics, and set the y2tics. We also set the y2label and xlabel. This latter one is empty, but we still need the space for it, so we set it to ' '. Beware the white space!

The next important step is the setting of the aspect ratio of our figure, which will be 0.6:1. Having done that, we introduce 3 labels: one for the xlabel, and two for the key. The placement of these labels is somewhat arbitrary, and depends on the particular terminal that you use. But only these three numbers and the coordinates of the key need any tweaking, really. At the very end, we plot our data. I plotted the second column twice, with an added random numbers, so that we can have two columns at each data point. Note that we plot the first column, too, but pass it to the xlabel command. By doing that, we can automate the labelling of the xtics, using the first column in our data file. Also note the specification of the titles: in the first two cases, the single quotes include a white space, while in the third one, there is nothing.

The impossible gnuplot graphs

2009-09-15T14:52:00.000-07:00

This post is going to be a short one. I would only like to announce that most of the material presented here is collected in a more orderly fashion on my web page .

Comments should still be posted here.
Cheers!

The slice of the pie II.

2009-09-13T09:14:00.000-07:00

Yesterday, I tried my hand at an exploded 3D pie chart. While we can use that method in particular cases, it is not "fool-proof": we might have to set the parameters by hand, should something cover something that it wasn't supposed to cover. Obviously, this situation is not ideal, even a bit disturbing, so I began to think what we could do differently.

The problem yesterday was that we used multiplot, and multiplot does not respect the depth ordering, i.e., we cannot let gnuplot determine the distance of the elements from the viewer, for it would miserably fail, simply because those elements are not parts of the same plot. Therefore, this whole affair could be much easier, if we were able to manage to put everything, or at least, the crucial pieces in one plot. Well, I have done the hard work, you have just got to read on to see the magic trick. I will start out with a particular example, but then I will show how this can be automated, so you just have to give the file name, and the rest will be done by the script. But, just to wet your appetite, here is the figure

Now, let us get down to business. We have six data points, which, for the sake of simplicity, are called

A=0.05*2*pi; B=0.3*2*pi; C=0.4*2*pi; D=0.1*2*pi; E=0.1*2*pi; FF=0.05*2*pi;

The last one is not really relevant, for the sum of these numbers is 2 pi(e), and FF is not used anywhere at all. If you recall, the reason for having to use multiplot was that when we put the pieces together, we changed the palette after each step, so the next piece, by default, had to be in a separate plot, and then we couldn't use depth ordering any more. Since we need different colours for the slices, the only way out of the loophole that I mentioned above is that we plot everything into a file, and then plot the file. When doing so, we need to specify the colour by a separate function, but that is really easy. After this introduction, let us see the script, which I will discuss further.

reset
b=0.5; a=0.5; r=1.0; s=0.1; m=1.5; eps=1e-4; N=6
A=0.05*2*pi; B=0.3*2*pi; C=0.4*2*pi; D=0.1*2*pi; E=0.1*2*pi; FF=0.05*2*pi

f(x,n) = (x>n?0.0:1.0)
F(x) = 1+f(x,A-eps)+f(x,A+B)+f(x,A+B+C)+f(x,A+B+C+D)+f(x,A+B+C+D+E)
at(y,x) = (x==0.0?0:(atan2(y,x)>0.0?F(atan2(y,x)):F(atan2(y,x)+2*pi)))

c(u,v,q)=cos(u)*r*v+q*cos(A/2); s(u,v,q)=sin(u)*r*v+q*sin(A/2)
C(u, q) = cos(u)*r+q*cos(A/2); S(u,q) = sin(u)*r+q*sin(A/2)
z = s+a; Z(v) = s+a*v

rg(x) = abs(cos(100*x/7))
gg(x) = abs(cos(100*x/11-pi/2))
bg(x) = abs(cos(100*x/13+pi/2))

set view 30, 20; set parametric
unset border; unset tics; unset key; unset colorbox; set ticslevel 0
set pm3d depthorder; set pal maxcolor N+1
set vr [0:1]; set xr [-1.5:1.5]; set yr [-1.5:1.5]; set zr [0:2]; set cbr [0:2*pi]

set multiplot
set iso 2, 2
set table 'pieslice2.dat'
set ur [A+eps:2*pi]
splot C(u,-eps), S(u,-eps), Z(v)
set ur [0+eps:A-eps]
splot C(u,b), S(u,b), Z(v)
set ur [0:1]
splot u*r, 0, Z(v), u*r*cos(A), u*r*sin(A), Z(v), \
u*r+b*cos(A/2), b*sin(A/2), Z(v), u*r*cos(A)+b*cos(A/2), u*r*sin(A)+b*sin(A/2), Z(v)
unset table

set iso 2, 80
set table 'pieslice3.dat'
set ur [A+eps:2*pi]
splot c(u,v,0-eps), s(u,v,0-eps), z
set iso 2, 2
set ur [0:A]
splot c(u,v,b), s(u,v,b), z
unset table
unset multiplot

set multiplot
set pal functions rg(gray)/m, gg(gray)/m, bg(gray)/m
sp 'pieslice2.dat' u 1:2:3:(at($2,$1)) w pm3d

set pal functions rg(gray), gg(gray), bg(gray)
splot 'pieslice3.dat' u 1:2:3:(at($2,$1)) w pm3d
unset multiplot

OK, so let us see what is happening here! In the first line, we define a couple of things that will determine the look-out of our figure. 'b' will be the shift of the cut-out, 'a' is the pie's height, 'r' is its radius, 'm' will determine the shade of the sides, with respect to the top, and 'eps' is just a small number whose role will become clear in a second. And finally, 'N' is the number of our data points. Then come the data points. After this comes the heavy part. Literally.

f(x,n) is a Heaviside function, which we use in the colour scheme. Remember that we want to colour the pie according to a function that depends on the azimuth angle. We will have a number of colours, and change the pm3d colour whenever the angle crosses the next value in our data set. If you look closely, our F(x) increases by one at each such point. (This is why we needed the Heaviside function.) That is, we could use F(angle) to colour our plot! However, there is a small glitch: since we are going to plot into a file first, we lose the information on the angle, and we will have to undo the polar-Cartesian transformation. This is why we need some distorted version of the atan2 functions. With the definition given here, we can make sure that it returns values in [0:2 PI], and there isn't a phase jump at (-1, 0).

Having defined these 3 functions, we define our shapes. These are identical (with the exception of the trivial third/second variable, 'q') to those shown yesterday, so I will not discuss them here. The next three lines will determine the colour scheme. If you are unhappy with the pie that you get, you should change these lines here.

The next 4 lines define some properties of the image. There are only two things to watch out for: one is that we fix the number of colours in our palette, namely, 7. The second is that we fix the range of the palette: this is done by the

set cbr [0:2*pi]

where cbr stands for cbrange.

By now, we have set everything, so we are ready to plot, first to files. We will have two files, pieslice2.dat and pieslice3.dat. The first one will hold the side of the pies, while the second is the top. The only reason for putting them in two separate files is that we want darker colours for the sides. These plots are trivial, and should speak for themselves. The only exception is the shift by 'eps'. We did this, so that the atan2 function is everywhere defined. Otherwise, we would have points at (x,y) = (0,0), where the atan2 function would return an undefined value, and that would just give us a hole in the plot. Also note that we have set the number of isosamples. We have only 2, where high resolution is not required, and set it to larger numbers only where needed. This improves on speed, but more importantly, reduces the size of the file considerably, if redirected to some vector graphics format, e.g., eps, or pdf. Believe me, this really makes a lot of difference!

At this point, we have all data in two files, so we can simply plot them. We plot them separately, because we want to colour them differently, so we have to change the palette. But these steps are really straighforward. The pie is there, if you want to put a background, labels, etc., you can do it here.

At the beginning, I mentioned that this whole business can be done automatically. If you look at the script above, you will notice that the variables and all related functions are defined at the beginning: N was the number of samples, then we had the data points, and finally, F(x) depends on the data. So, if we have a script that prints all these into a file, we can load that file there, and the rest is unchanged. The improved version would look like this:

reset
b=0.5; a=0.5; r=1.0; s=0.1; m=1.5; eps=1e-4

f(x,n) = (x>n?0.0:1.0)
load 'pie_l.gnu'
at(y,x) = (x==0.0?0:(atan2(y,x)>0.0?F(atan2(y,x)):F(atan2(y,x)+2*pi)))

c(u,v,q)=cos(u)*r*v+q*cos(A/2); s(u,v,q)=sin(u)*r*v+q*sin(A/2)
C(u, q) = cos(u)*r+q*cos(A/2); S(u,q) = sin(u)*r+q*sin(A/2)
z = s+a; Z(v) = s+a*v

rg(x) = abs(cos(100*x/7))
gg(x) = abs(cos(100*x/11-pi/2))
bg(x) = abs(cos(100*x/13+pi/2))

set view 30, 20; set parametric
unset border; unset tics; unset key; unset colorbox; set ticslevel 0
set pm3d depthorder; set pal maxcolor N+1
set vr [0:1]; set xr [-1.5:1.5]; set yr [-1.5:1.5]; set zr [0:2]; set cbr [0:2*pi]

set multiplot
set iso 2, 2
set table 'pieslice2.dat'
set ur [A+eps:2*pi]
splot C(u,-eps), S(u,-eps), Z(v)
set ur [0+eps:A-eps]
splot C(u,b), S(u,b), Z(v)
set ur [0:1]
splot u*r, 0, Z(v), u*r*cos(A), u*r*sin(A), Z(v), \
u*r+b*cos(A/2), b*sin(A/2), Z(v), u*r*cos(A)+b*cos(A/2), u*r*sin(A)+b*sin(A/2), Z(v)
unset table

set iso 2, 80
set table 'pieslice3.dat'
set ur [A+eps:2*pi]
splot c(u,v,0-eps), s(u,v,0-eps), z
set iso 2, 2
set ur [0:A]
splot c(u,v,b), s(u,v,b), z
unset table
unset multiplot

set multiplot
set pal functions rg(gray)/m, gg(gray)/m, bg(gray)/m
sp 'pieslice2.dat' u 1:2:3:(at($2,$1)) w pm3d

set pal functions rg(gray), gg(gray), bg(gray)
splot 'pieslice3.dat' u 1:2:3:(at($2,$1)) w pm3d
unset multiplot

where 'pie_l.gnu' contains the following two lines

N=6
F(x) = 1+f(x,0.314)+f(x,2.199)+f(x,4.712)+f(x,5.34)+f(x,5.969)

and nothing else. At this point, we can either let a script write these numbers to 'pie_l.gnu', or just print it to the standard output, and re-direct that output to gnuplot as

load '< somescript somedata.dat'

The following gawk script should do the job, provided that you want to process the first column in your data file. The script also normalises, so arbitrary numbers can be used.

#!/bin/bash

gawk 'BEGIN {sum=0.0; i=0}
    {       a[i] = sum+$1
            sum += $1
            i++
    }
    END {   printf "N=%d\n", i
            printf "F(x) = 1"
            for(j=0;j<i-1;j++) printf "+f(x,%f)", 6.28318530717959*a[j]/sum
            printf "\n"
    }' $1

This was a long post, but I hope that it is more or less clear what we are doing, and after all, it is not hard at all.
Cheers!

The slice of the pie

2009-09-12T11:48:00.000-07:00

I have already shown various ways of making a pie chart in gnuplot. We will brush off one of the methods, and see how we can move one of the slices a bit. At the end of the day, we will have the following graph

For the sake of simplicity, I will demonstrate the method for four data points, which we will call

A=0.2*2*pi; B=0.3*2*pi; C=0.4*2*pi; D=0.1*2*pi;

In this particular case, they add up to 2 pi(e);) Once we see what and how to do, we can drop this condition, and use any data. But first, let us see the script

reset
b=0.5; a=0.5; r=1.0; s=0.1; m=1.5
c(u,v)=cos(u)*r*v+BB*cos(A/2); s(u,v)=sin(u)*r*v+BB*sin(A/2)
C(u) = cos(u)*r+BB*cos(A/2); S(u) = sin(u)*r+BB*sin(A/2)
z = s+a; Z(v) = s+a*v
A=0.2*2*pi; B=0.3*2*pi; C=0.4*2*pi; D=0.1*2*pi;
set view 30, 20; set parametric
unset border; unset tics; unset key; unset colorbox
set ticslevel 0
set vr [0:1]; set xr [-1.8:1.8]; set yr [-1.8:1.8]; set zr [0:2]
set multiplot
set ur [0:1]; set pal mo RGB func 0.8, 0.8, 0.8
splot 3.6*u-1.8, 3.6*v-1.8, 0 w pm3d

BB=b
set ur [0:r]; set pal mo RGB func 1/m, 0, 0
splot BB*cos(A/2)+cos(A)*u, BB*sin(A/2)+sin(A)*u, Z(v) w pm3d
splot BB*cos(A/2)+u, BB*sin(A/2), Z(v) w pm3d

set ur [0:A]; set pal mo RGB func 1/m, 0, 0
splot C(u), S(u), Z(v) w pm3d
set ur [A:A+B]; set pal mo RGB func 0, 1/m, 0; BB=0.0; rep
set ur [A+B:A+B+C]; set pal mo RGB func 0, 0, 1/m; rep
set ur [A+B+C:A+B+C+D]; set pal mo RGB func 1/m, 1/m, 0; rep
set ur [0:r]; set pal mo RGB func 0, 1/m, 0
splot cos(A)*u, sin(A)*u, v*a+s w pm3d

BB=b
set ur [0:A]; set pal mo RGB func 1, 0, 0
splot c(u,v), s(u,v), z w pm3d
set ur [A:A+B]; set pal mo RGB func 0, 1, 0; BB=0.0; rep
set ur [A+B:A+B+C]; set pal mo RGB func 0, 0, 1; rep
set ur [A+B+C:A+B+C+D]; set pal mo RGB func 1, 1, 0; rep
unset multiplot

At the beginning, we define several plot-related constants, such as, 'b', which gives the displacement of the red slice, a, which is the height of the pie, r, which is the radius or the arcs, s, which determines by how much the pie is lifted from ground, and m, which we will use in the colouring: the sides of the slices will be m times darker, than the top, so as to give the impression that the top is lit, while the sides are in the shadow.

Then we define 6 functions, which are just helpers, so that the rest of the code will look tidier. These are basically the parametrisation of the sides and the tops. Having done this, we define A, B, C, and D, our data points, then we set a couple of things on the figure, and in the multiplot, plot the "ground".

Once we have the ground on which to build, we plot the sides of the pie. Note that since we have two cuts, we would, in principle, need 4 new surfaces. However, one can never be seen, so we can just drop that altogether. Then we plot the cylinder, each segment in a different colour. Since those plots are actually the same, except for the range that we re-set every time, we can simply replot everything. This makes the code a bit cleaner, I believe. Also note that the colours are divided by m, which we defined at the beginning. Finally, we put on the cap, in a brighter colour. Adding the labels, if needed, should be quite straightforward. You can look up the details in my other pie chart posts.

The only trick that we used here was that we shifted the cut-out slice by such an amount that only one surface was partially covered by the rest of the pie. In general, it is quite hard to determine which surface will cover which, if we let the cut-out closer, so I cheated here a bit.

Bubble graphs in gnuplot, once more

2009-09-03T02:24:00.001-07:00

We have already discussed two ways of making a bubble plot. Now we will see a third one, which is probably the best of these, for it doesn't have the overhead of a) hacking the postscript file, b) plotting something with pm3d. This is so simple! Here is my graph

and here is the code

reset
set sample 30
set table 'bubble.dat'
p [0:5] sin(x*3.0)*exp(-x)+rand(0)/10.0
unset table

xs=-0.007
ys=0.002
unset key
set object 1 rect from screen 0,0 to screen 1,1 behind fillcolor rgb "#ddddff"
set border back
set xlabel 'Time [s]'
set ylabel 'Position [m]'
p [-0.1:5] 'bubble.dat' u 1:2 w p ps 3 pt 7 lc rgb "#ff0000", \
'' u ($1+xs):($2+ys) w p ps 2.6 pt 7 lc rgb "#ff2222", \
'' u ($1+xs):($2+ys) w p ps 2.2 pt 7 lc rgb "#ff4444", \
'' u ($1+xs):($2+ys) w p ps 1.8 pt 7 lc rgb "#ff6666", \
'' u ($1+xs):($2+ys) w p ps 1.4 pt 7 lc rgb "#ff8888", \
'' u ($1+xs):($2+ys) w p ps 1.0 pt 7 lc rgb "#ffaaaa", \
'' u ($1+xs):($2+ys) w p ps 0.6 pt 7 lc rgb "#ffcccc", \
'' u ($1+xs):($2+ys) w p ps 0.2 pt 7 lc rgb "#ffeeee"

The red parts are really irrelevant: the first one just produces some dummy data, we have seen this so many times that it is becoming boring. Then we define two variables (blue) 'xs' and 'ys', and draw a rectangle, just to have some background. In the last 8 lines we plot our data 8 times, each time in a different colour (the colour successively becomes whiter) and in different size (the points successively become smaller). At the same time, the centre of the circles is shifted by 'xs' and 'ys', so that we give the impression that the spheres are lit from the top left corner. By modifying these values, you can move your virtual light source.

The are only three caveats here: one is that since we apply a shift to the data points, we have got to make sure that our xrange and yrange supports the shifter data. This is why, while we have data in [0:5], our xrange is actually [-0.1:5]. The yrange does not require special attention in this case.

The second caveat is that it is perfect for a raster image, it will look quite nasty on a vector image. If you plan on printing the graph through postscript or pdf, you will have to define more levels for the transition from red to white. But the idea is the same, you will simply have more plots.

And third, when you define your 'xs' and 'ys', you have to make sure that the successive plots fall on the first one, otherwise your bubbles will have some funny shape.

Filled curves in gnuplot

2009-08-31T12:19:00.000-07:00

Gnuplot has a plotting option called filledcurve. I do not want to dwell on the specifics of this, because I think that most about it can be learned by looking at the demo page. What I would like to discuss instead is how this can be used to turn our graph into a sleek figure. This is not something that you would put in a paper publication, on the other hand, on a poster, a web page or similar, it looks much better, than what you would get using the default. I just show the image, and we will discuss the script later.

reset
f(x) = x+20
g(x) = x*x*x-2.0*x*x+x+10-x*x*x*x*0.01
set sample 101
set parametric
set trange [0:4]
set table 'filled.tab'
plot g(t)+5.0*(0.5-rand(0)), f(t)+5*(0.5-rand(0))
unset table
unset parametric
set iso 2, 100
set xrange [0:100]; set yrange [0:50]
set table 'filledbg.tab'
splot y*0.5
unset table

set key left
unset colorbox
set border lc rgb "white"
set grid front lc rgb "#dddddd"
set xlabel 'Time [a.u.]' tc rgb "white"
set ylabel 'Amplitude [a.u.]' tc rgb "white"
set palette defined (0 0.2 0.2 0.2, 1 0.6 0.6 0.6)
set object 1 rect from screen 0, 0 to screen 1, 1 behind fc rgb "black" fillstyle solid 1.0
p 'filledbg.tab' w image t '', \
'filled.tab' u 0:1:2 w filledcurve above lt 1 lc rgb "#467F1E" t '', \
'' u 0:1:2 w filledcurve below lt 1 lc rgb "#802020" t '', \
'' u 0:1 w l lt 8 lw 1 t 'First', '' u 0:2 w l lt 9 lw 1 t 'Second'

The first part of the code is nothing but generating some dummy data. If you have something to plot, this is unnecessary, and the actual plot takes only 12 lines altogether. However, if you want to have a gradient in the background, you will have to generate 'filledbg.tab'. This is the file that we will plot as an image, and this is what will give the gradient. I would also point out here that in order to generate 'filledbg.tab', we need the range of the actual plot. If you don't know it beforehand, you can use a trick to determine this. I have discussed this at great length in a few of my previous posts, so you could look it up there.

Now, sometimes you do need a background for the plot, and this is not just some silly whim, but you might have to match the background of your poster, or your transparency. In order to achieve this, we will draw a rectangle behind the graph, and fill it with the desired colour. What we should keep in mind, however, is that without taking extra measures, the axis labels will have the default colour of black, and this will not look very good on a black background. Therefore, you should specify the colour in the declaration of the labels.
In the plot, first we draw the gradient. Without defining the grid in the front, the grid would be obscured, therefore, we push it to the forefront. The real curves are drawn four times: first we draw the filled curves, simply using three columns of our data file. (One column can be dropped, in which case an axis or a function could be used to define the boundaries.) If the two curves cross each other, we could have different colours, depending on the relative size of the functions that we plot. This is why we have to plot our curves for the second time: we have different colour for when the first curve is above and below the second one. The last two plot are just there to draw the actual curves in different colours.

The plot thickens - shadow to a curve

2009-08-30T06:16:00.000-07:00

This is Sunday afternoon, so I will discuss an easy plot, namely, how we can give the impression that a 2D plot is in 3D, simply by adding a slight shadow to a curve. And the plot thickens, because the curve itself is thick.:) Here is the figure that we will have

and here is the script that produced the image

reset
h(x) = x*x*x-5*x*x+x+20
a = 0.05
b-0.5
unset colorbox

set iso 2, 2
set table 'thickbg.dat'
sp [0:5] [0:30] y
unset table

set sample 20
set table 'thick.dat'
plot [0:5] h(x)+5*(0.5-rand(0))
unset table

set xlabel '{/Helvetica=14 Time}
set ylabel '{/Helvetica=14 Price}'
set tics font "Helvetica, 14" nomirror
set grid front
unset grid
set palette defined (0 0.95 0.95 0.95, 1 0.95 0.95 0.95)
plot [0:5] [0:30] 'thickbg.dat' w ima t '',\
'thick.dat' u ($1+a):($2+b) w l lc rgb "#eeeeee" lw 19 t '',\
'' u ($1+a):($2+b) w l lc rgb "#dddddd" lw 17 t '',\
'' u ($1+a):($2+b) w l lc rgb "#cccccc" lw 15 t '',\
'' u ($1+a):($2+b) w l lc rgb "#bbbbbb" lw 13 t '',\
'' u ($1+a):($2+b) w l lc rgb "#aaaaaa" lw 11 t '',\
'' u ($1+a):($2+b) w l lc rgb "#999999" lw 9 t '',\
'' u 1:2 w l lc rgb "#ff0000" lw 10 t  ''

The whole (and very simple) idea is that we will add a replica of the curve that we want to plot to the figure, in some shade of gray, and shifted a bit, so as to give the impression that this curve is the shadow of the other one.

First, we define three functions. h(x) will generate some data (if you have a data file to plot, then obviously, you can skip this, and the set table... plot h(x) ... unset table section), and 'a' and 'b' define the shift in the horizontal and vertical directions, respectively. We also add a background to the figure, in gray. In order to do so, we splot 'y' (or anything else, for that matter) to a file, and in our real plot, we plot this as an image, with the defined colour palette.

There is some subtlety to the lines

set grid front
unset grid

At first, this appears silly to set something and then unset it immediately afterwards, but there is a very good reason for that. Odd as it might sound, the placement of the tics is defined by set grid, i.e., if we want visible tics on the background, we have to push them to the front, otherwise, the background of our figure would cover them. We do this by setting grid front. However, we do not actually want any grid lines, so we unset them. The effect on the tics remains, however.
The last step before plotting our curve is to plot the shadow. We do it in 6 steps, and in each step we reduce line width and the whiteness of the colour. If you look at the colour definitions, consecutive lines are in darker shades of gray. The reason for this is that in this way the edge of the shadow will not be sharp. However, if you are satisfied with a sharp shadow, you can skip 5 out of these 6 steps, and draw the shadow only once. In fact, this does not look that bad at all, as you can see in the figure below.

Well, this is all for today.

Changing the aspect ratio of axes

2009-08-21T11:58:00.000-07:00

A couple of days ago, I saw a question on the gnuplot discussion board about how to change the ratio of the three axes in a 3D plot independently. In 2D, this is not a problem, for one can use the command

set size xsize, ysize

where xsize and ysize are relative to the canvas size. If, e.g., you type

set size 1, 0.5

then the figure will be just as wide as it was originally, but the height will be reduced by a factor of two. Now, the problem is that in 3D the above-mentioned instruction will lead to a figure whose height (i.e., the z coordinate axis) is reduced by a factor of two, but the width is unchanged. However, the width is determined by x and y, therefore, the x-y axes will be linked in this way. This results in the pitiable situation that we can't have an x axis, which is twice as long as the y axis. There's got to be a way around this problem! Sure enough, there is, and this is what we will set out to solve now. This will require some hand-work, but it is fairly simple otherwise.

The original problem on the discussion board was to set the xrange as [-6:6], the yrange as [-3:3], and the zrange as [-2:2], while keeping the units on each axis equal, i.e., have aspect ratios 6:3:2. What we can do is the following: we can easily change the z:x ratio, simply setting the size. Then, instead of plotting over the whole y range, we will only plot over one part of it, and yes, you've guessed correctly, we will use the ternary operator for that. The only trick is that we specify a yrange which is actually larger than the one that we want to plot. So, here is our script, and we will discuss it below

reset
set iso 100,100
set xrange [-6:6]
set yrange [-3:9] # We should have [-3:3], but do not rush!
set zrange [-2:2]
set ticslevel 0
unset key
unset colorbox
set xtics nomirror
set ytics -3,2,3
set ztics -2,2,2
set xtics -6,2,6
set border 1+16

set xlabel 'x-label'
set zlabel 'z-label' offset graph 0.05,0,-0.1
set ylabel 'y-label' offset graph -0.05,-0.25,0
f(x) = (x<3?x:1/0)
set size 1,0.7
set arrow 1 from -6,-3,-2 to -6,3,-2 nohead
set arrow 2 from 6,-3,-2 to 6,3,-2 nohead
set arrow 3 from -6,3,-2 to 6,3,-2 nohead

set palette rgbformulae 33,13,10
splot sin(x)*cos(f(y)) with pm3d

In the beginning, as usual, we set up the figure. The only thing to watch out is the unusual yrange: we specify more than we actually will use. On the border, we draw only the x, and z axes, but not the x2, y, and y2 axes. We will do that by hand, where we draw 3 headless arrows. Before we plot anything, we define a function which we will use to restrict the plotted yrange to [-3:3]. In the plotting function we pass not y, but f(y), so that nothing will be plotted beyond y=3. Having done the heavy work, we plot the function, and get the following image

Now, a couple of comments are in order here. One is that, since we set the y axis by hand, we have got to set the ylabel by hand, too. This is why we specified the offsets for the labels. (For the zlabel it wouldn't have been necessary, it just made the figure look a bit tidier.) The second remark is that it really depends on the view what looks the same on the 3 axes: that is why there is no generally applicable number can be given in the set size command. You have to try and determine it on a case-by-case basis.

Plotting the recession - some basic data analysis with gnuplot

2009-08-09T12:44:00.000-07:00

Today, we will see how we can do some very rudimentary data analysis, using gnuplot. Of course, if you want to analyse your data, then you can either use a dedicated script, written in a far better fit language, e.g., in Perl or in gawk, or you can do the whole thing in a tabular environment, using Excel or Origin, or Grace, or Labplot, or I do not know what. And I have also read somewhere that gnuplot conforms with the Unix philosophy, "one task - one tool", and we should not do any data processing with it. Now, there are at least two problems that I see with this approach: one is that some data analysis are really simple, and it is not clear, why one should use a separate tool for that. Furthermore, it will always cause some trouble (at least, in my experience) when the data analysis is not done from gnuplot, for one can never be sure how and on which version of the analysed data file the figure was based. And pipes can only be used on Linux/Unix-type systems, thus, the platform independence of gnuplot is lost. But more importantly, and if we are at the "philosophy" anyway, it is completely unclear (at least, to me) what constitutes plotting: just placing the curves, points and dots on an empty canvas, and writing a couple of labels for the axes, or this and the choice of the curves, points and dots. The problem is that one can in most cases decide what to do with erroneous points and the like only when it is already plotted and shown. So, if that is the case, why should we sequester data analysis from the plotting? Are they not the two complementing parts of the very same procedure?

So, these are my reasons for doing data processing in gnuplot. Of course, in most cases, I do that using an external script, and piping it from the gnuplot prompt. But that would deprive windows users from the pleasures of doing something fancier:)

By the way, if you look at some of the demo scripts of gnuplot 4.3, you will realise that the developers of gnuplot have already parted with the idea of "one task - one tool".

The kind of data processing we will do today requires a "recursive" access to data points, in particular, I will demonstrate how we can calculate the derivative and the integral of a series of data points. Calculating the derivative and integral of an analytic function is trivial, even if you do not have the analytical answer: for the derivative, you can simply plot f(x+dx)-f(x), where dx is just some pre-defined number, and for the integral you can apply a recursive scheme, detailed on the gnuplot demo web site. But they work on known functions, and we would like to calculate something based on "measured" data. For this, I will take the weekly oil price from this year, which you can download from the link. Let us suppose that we plot the oil price, but we also want to highlight the time intervals where the price is dropping. This requires us to calculate the derivative at each point (In this particular case, we could get away just by comparing the value and the next, but we have got to get access to the values anyway, and we would get it by fitting a linear function, whose linear coefficient gives the derivative. So, for better of worse, I will just stick with the derivative.)

For a starter, let us see our two scripts! As I have already pointed out, we will have to walk through the values, one by one, so we will need a 'for' loop, which means two scripts: a "master", and one that we call repeatedly. I have discussed this method at great length in my recent posts.

reset
unset key
A=0; sa=1; sb=0; integ=0.0
f(x) = a*x+b
g(x) = (x>0?GPVAL_Y_MIN:GPVAL_Y_MAX)
set xlabel 'Time [weeks]'
set ylabel 'Crude oil price [$]'
p 'oil.dat' u 0:8

xmax = GPVAL_DATA_X_MAX

set print 'oild.dat' append
l 'derint_r.gnu'
set print

plot 'oild.dat' u 1:(g($2)) w filledcurve lc rgb "#eeeeee", 'oil.dat' u 0:8 w l

and our 'for' loop, that we call 'derint_r.gnu'

a=1.0; b=40.0
fit [A:A+1] f(x) 'oil.dat' u 0:8 via a, b
integ = integ+(a*A+b+a*(A+1)+b)/2.0
if(sa*a<0) print A, sa, sb, integ; print A, a, b, integ; sa=a; sb=b
if(sa*a>0) print A, a, b, integ; sa=a; sb=b
A=A+1
if(A<xmax) reread

The first part of the main script should be clear: we set up the style of our figure, and do a dummy plot to learn something about our data file, 'oil.dat'. f(x) is our linear function, while g(x) is just a helper, whose role will become evident a bit later. Having set up the plot, we call 'derint_r.gnu' 'xmax' times. So, what is in 'derint_r.gnu'? We simply fit our linear function over a range that contains only two points, and we step through the whole file by moving those two points by one in each call to the loop. Then we calculate the integral by adding the area of the next trapezium, and then depending on whether the sign of the derivative, 'a', is different to the previous one, we print out the results once (if the signs are identical), or twice (if they are different). This distinction is necessary because we want to plot rectangles where the derivative is negative. After we are done with the data processing, we close our new file by setting the print, and plot the results, to get the following figure

Note that in the plotting, we called the function g(x), which takes on the value of the minimum of the yrange, if the argument is negative (the price is dropping), while returns the maximum of the yrange otherwise. Of course, there are many things that we could do to improve the image, but I wanted only to demonstrate how we can do the data processing part. You can dress up the image as you like, and for that you can also draw on many a post of mine.

Return of the wall chart - entirely in gnuplot

2009-08-08T10:59:00.000-07:00

On the 7th of June, I showed how one can draw a wall chart in gnuplot. However, that method relied on the knowledge of the particular function we wanted to plot, therefore, would not work for data files. We have since learnt a thing or two, and it is high time to return to that type of plot, and see whether we could use something that we have picked up on the way. Yes, you are right, the answer lies in the affirmative, and the situation is even far better: all we have to do is to take the scripts from yesterday, and change two lines. It could not be simpler!

Before we turn to the complete chart, first we should see how we can plot just one single plane, taken the values from the file. For the sake of simplicity, I will use 'ribbon.dat', given in the post yesterday. Also, if you haven't read that post, you should now, for I will simply change two lines, but I will not go through the whole code again.

If you recall, our data file, after some processing in gnuplot, looked like this (well, at least, the first column)

#Surface 0 of 1 surfaces

#IsoCurve 0, 3 points
#x y z type
0  13  0.901743 i
1  13  0.901743 i
2  13  0 i

#IsoCurve 1, 3 points
#x y z type
0  12  0.860648 i
1  12  0.860648 i
2  12  0 i
...

This is what we called 'rib2.dat'. We used the every keyword of splot to single out the first (0) and second (1) line in gnuplot. Now, what would happen, if we plotted the second and third lines? We would get a single plane, standing vertically, and having a "shape" given by the 'z' values, as in this figure:

Just for the record, our plot is produced by these three lines in 'ribbon_r.gnu'

r=rand(0); g=rand(0); b=rand(0)
set palette model RGB defined (0 r/cm g/cm b/cm, 1 r/cm g/cm b/cm)
splot 'rib2.dat' every ::1::2 u (A-1):2:3 w pm3d

i.e., we restrict the 'x' values to 'A-1', and plot the second and third columns of the second and third records in each block. (Remember that the numbering of the records begins with 0!)
The figure is almost OK, with the tiny glitch that the negative values are plotted "downwards", i.e., the reference line is 0. We can fix this with the helper function

zr(x) = (x==0?zmin:x)

which returns the minimum of the zrange, if the 'z' value is smaller, than 0, otherwise it returns 'z'. So, by modifying the plot as

r=rand(0); g=rand(0); b=rand(0)
set palette model RGB defined (0 r/cm g/cm b/cm, 1 r/cm g/cm b/cm)
splot 'rib2.dat' every ::1::2 u (A-1):2:(zr($3)) w pm3d

we get

Remember that, since we use a random number to determine the RGB values, the colour scheme of the plot will be different after each call of the script, and this is why the plot looks a bit different to the previous one. Nevertheless, the point is not this, but that the plots go up from the bottom of the graph to the corresponding 'z' values. If you are satisfied with this, you can stop here, there is nothing more to it. However, we could turn what we have into a "real" wall chart by adding a "roof" to the walls. All it takes is four extra lines in our 'for' loop, which will now look like this (this is wall_r.gnu for reference)

A=A+1
set table 'rib.dat'
plot filename u A:A
unset table

set table 'rib2.dat'
splot 'rib.dat' mat
unset table

r=rand(0); g=rand(0); b=rand(0)
set palette model RGB defined (0 r/cm/cm g/cm/cm b/cm/cm, 1 r/cm/cm g/cm/cm b/cm/cm)
splot 'rib2.dat' every ::1::2 u (A-1):2:(zr($3)):3 w pm3d

set palette model RGB defined (0 r g b, 1 r g b)
splot 'rib2.dat' every ::0::1 u (A-1+B*$1):2:3 w pm3d

set palette model RGB defined (0 r/cm g/cm b/cm, 1 r/cm g/cm b/cm)
splot 'rib2.dat' every ::1::2 u (A-1+B):2:(zr($3)):3 w pm3d

if(A<C) reread

Note that we plot the same file, 'rib2.dat', thrice, each time with a different shade of the same colour: first with the darkest, and this is the plane farthest from the viewer, then the roof of the walls, i.e., the ribbons, with the lightest, and finally, the front plane of the walls with a medium shade. The main script reads as follows

reset
filename="ribbon.dat"
mult=1.2; cm=1.4
A=0; B=0.3;
set xlabel 'x axis [a.u.]'; set ylabel 'y axis [a.u.]';
set ticslevel 0
unset key; unset colorbox
set tics out nomirror
set border 1+2+4+8+16+32+64+256+512 back
mini(x)=(x<0?x*mult:x/mult)>0?x*mult:x/mult)
zr(x) = (x==0?zmin:x)

set isosample 100, 3
set xrange [0:1]; set yrange [0:1]
set table 'bg.dat'
splot x
unset table; set xrange [*:*]; set yrange [*:*]

splot filename mat
C=GPVAL_DATA_X_MAX+1
ymax = GPVAL_DATA_Y_MAX+1
zmin = mini(GPVAL_DATA_Z_MIN)
zmax = maxi(GPVAL_DATA_Z_MAX)
set xrange [-0.5:C-0.5]
set yrange [0:ymax]
set zrange [zmin:zmax]
set cbrange [0:1]

set multiplot
set palette model RGB defined (0 0.3 0.3 0.85, 1 0.9 0.9 0.95)
splot 'bg.dat' u ($1*C-0.5):($2*ymax):(zmin):3 w pm3d, \
'' u ($1*C-0.5):(ymax):(1.9999*$2*(zmax-zmin)+zmin):3 w pm3d, \
'' u (-0.5):($1*ymax):(1.9999*$2*(zmax-zmin)+zmin):(0) w pm3d
set cbrange [zmin:zmax]

unset border; unset tics; unset xlabel; unset ylabel; unset zlabel

l 'wall_r.gnu'

unset multiplot

This is identical to 'ribbon.gnu' from yesterday, with the exception of the definition of the helper function 'zr(x)'. The resulting figure looks like this one

There are two more things that we could do with this graph: One is that we can add a grid to the surfaces by adding these three lines to our 'for' loop

set style line 1 lt -1 lw 0.5 ps 0
set pm3d implicit at s
set pm3d depthorder hidden3d 1

(Or, we could add this to our main script, but only after plotting the background, otherwise, that will be gridded, too.) By doing so, we get the following figure

The second is that instead of drawing the back plane of the wall, we could add a front plane, thereby "closing" the volume. We do this by replacing the first plot command in our 'for' loop by

splot 'rib2.dat' every ::0:(ymax-2):1:(ymax-1) u (A-1+B*$1):(0):(zz($2,$3)) w pm3d

and adding the helper function

zz(x,y) = (x==1?zmin:y)

to the main script. What the plotting command does is to take only the last two blocks, and choose only the first and second records. The last two blocks will have 1 and 0 in the second column, so we use that and the helper function to determine which 'z' value we want to pick to draw our rectangle in the front.
Of course, you should plot the front plane after you plotted the ribbon on the top, otherwise, parts of it might be covered. Once you do that, you get the following image

I believe, there is not too much else that we could do here, but we have gone a long way, and produced a decent-looking wall chart! I should also point out that, if you want to use the second plot, i.e., the wall chart without the roofs, you can still add a grid. All you have to do is to add those three lines that we discussed in connection with the grid on the complete wall chart, and you can, naturally, add a grid to the ribbon plot that we discussed yesterday.

Ribbon charts - without the gawk script

2009-08-07T09:00:00.000-07:00

Yesterday, we saw how we can produce a ribbon chart from an arbitrary data file, and with the help of an external script. You can still make it work under windows, although you might have to call the script outside gnuplot. Probably, it would be better to do everything in gnuplot, wouldn't it? So, can we do something to alleviate the situation. But, of course, we can, we just have got to use those little grey cells! If you are patient and read further, I will explain how you can draw this graph completely automatically! (Uhm, this is rather similar to that from yesterday... Never mind, the point is that you haven't got to do any work!)

If you recall, the key to plotting the ribbons was to convert the original data file (see the beginning of my previous post) into something like this:

0 0 0.857354
1 0 0.857354

0 1 0.551386
1 1 0.551386

0 2 0.243002
1 2 0.243002

0 3 -0.24
1 3 -0.24

0 4 -0.422674
1 4 -0.422674
...

i.e., we need two rows of three columns, separated by a linefeed. We will need some data processing for this, but we can manage. First, we will duplicate all columns, simply plotting into a file. By issuing

set table 'rib1.dat'
plot 'ribbon.dat' u 1:1
unset table

we get the following file, 'rib1.dat'

#Curve 0 of 1, 14 points
#x y type
0.857354  0.857354  i
0.551386  0.551386  i
0.243002  0.243002  i
-0.24 -0.24  i
-0.422674 -0.422674  i
-0.513352 -0.513352  i
-0.660823 -0.660823  i
-0.0962106 -0.0962106  i
-0.0152513 -0.0152513  i
0.489453  0.489453  i
0.876023  0.876023  i
0.945061  0.945061  i
0.860648  0.860648  i
0.901743  0.901743  i

This is almost what we want, two columns of the same numbers, namely those from the first column of 'ribbon.dat'. The only difficulty is that it contains a third column, which in this case is just a set of i, but if we were to plot this file as a matrix, at x=2 all values would be equal to zero, so it would be a distorted ribbon. How can we get rid of that letter? Well, we will turn the into a form similar to that at the beginning of this post, simply by

set table 'rib2.dat'
splot 'rib.dat' mat
unset table

where now 'rib2.dat' looks like this

#Surface 0 of 1 surfaces

#IsoCurve 0, 3 points
#x y z type
0  13  0.901743 i
1  13  0.901743 i
2  13  0 i

#IsoCurve 1, 3 points
#x y z type
0  12  0.860648 i
1  12  0.860648 i
2  12  0 i
...

which is very similar to the desired format, for now we can simply call an splot, without the matrix specifier. But we still have an erroneous value in every third line. We will remove that using the 'every' specifier of plot. If you haven't done it yet, you should definitely issue the command

?every

You will learn a lot of interesting things. What we want to do here is to plot only the first end second record in each block. In that case, it does not matter what is in the third record, because it will not be processed by gnuplot. So, we take our 'rib2.dat' file, and plot it as

splot 'rib2.dat' every ::0::1 u 1:2:3 w pm3d

which returns the first (0) and second (1) column in each block. Thus, with these three easy steps, we created a ribbon, which happens to represent the first column in our original data file. All that is left is to walk through the columns in 'ribbon.dat', and do the same thing again and again. Note that we specified the columns that we are going to use as 1:2:3. We will need this later, when we shift the next ribbon to a new position.

Then, let us suppose that we are a bit lazy, and we do not fancy following the prescription above. We can then resort to the recipe outlined a couple of days ago in connection with the pie charts and the bar graphs. If you haven't read those posts, you should now, for I am going to use everything mentioned in them.

Since the three steps above are repetitious, we will put them in a 'for' loop. As I discussed it earlier, the way to do this is to write the instructions in a separate file, and 'reread' them as many times as necessary. I.e., we will have a file, 'ribbon_r.gnu', which contains the following lines

A=A+1
set table 'rib.dat'
plot 'ribbon.dat' u 1:1
unset table

set table 'rib2.dat'
splot 'rib.dat' mat
unset table

r=rand(0); g=rand(0); b=rand(0)
set palette model RGB defined (0 r/cm g/cm b/cm, 1 r g b)
splot 'rib2.dat' every ::0::1 u (A-1+B*$1):2:3 w pm3d
if(A<C) reread

So, this is our 'for' loop, and this is where we call it, 'ribbon.gnu'

reset
filename="ribbon.dat"
mult=1.2; cm=2.5
A=0; B=0.3;
set xlabel 'x axis [a.u.]'; set ylabel 'y axis [a.u.]';
set ticslevel 0
unset key; unset colorbox
set tics out nomirror
set border 1+2+4+8+16+32+64+256+512 back
mini(x)=(x<0?x*mult:x/mult)>0?x*mult:x/mult)

set isosample 100, 3
set xrange [0:1]; set yrange [0:1]
set table 'bg.dat'
splot x
unset table; set xrange [*:*]; set yrange [*:*]

splot filename mat
C=GPVAL_DATA_X_MAX+1
ymax = GPVAL_DATA_Y_MAX+1
zmin = mini(GPVAL_DATA_Z_MIN)
zmax = maxi(GPVAL_DATA_Z_MAX)
set xrange [-0.5:C-0.5]
set yrange [0:ymax]
set zrange [zmin:zmax]
set cbrange [0:1]

set multiplot
set palette model RGB defined (0 0.3 0.3 0.85, 1 0.9 0.9 0.95)
splot 'bg.dat' u ($1*C-0.5):($2*ymax):(zmin):3 w pm3d, \
'' u ($1*C-0.5):(ymax):(1.9999*$2*(zmax-zmin)+zmin):3 w pm3d, \
'' u (-0.5):($1*ymax):(1.9999*$2*(zmax-zmin)+zmin):(0) w pm3d
set cbrange [zmin:zmax]

unset border; unset tics; unset xlabel; unset ylabel; unset zlabel
l 'ribbon_r.gnu'
unset multiplot

This looks formidable at first, but we can easily decipher all. First, in our 'for' loop, in 'ribbon_r.gnu', we simply wrote down the three steps required to created one ribbon, we re-set the colour palette, and call the loop as many times as many columns there are (C).

In the main script, in 'ribbon.gnu', we first set up the main properties of the figure, and we define the width (B) of the ribbons, the multiplier (mult), which determines by how much the figure is bigger, than the minimum and maximum of the data points in 'ribbon.dat', and we also define 'cm', which will determine how "deep" the colour modulation of the ribbons is. By setting this to a number close to 1, we get a more or less homogeneous colour, while setting it to something large, we get a wildly modulated colour palette for the ribbons. Note the two helper function, 'mini' and 'maxi', which are used to set the zrange, and which make use of the value of 'mult'. The next couple of lines produces the background, which you can skip, if you do not want a very fancy graph. Then we make a dummy plot, just to determine the various ranges. You can read about this in my last but one post. Immediately after the multiplot, we plot 'bg.dat', i.e., the background of the figure. After this we call 'ribbon_r.gnu', and finally, escape from the multiplot.

A couple of notes: in the 'for' loop, we generated the colours by picking three random numbers. This means that your plot will have different colours after successive calls of 'ribbon.gnu'. You can change this, if you want to have some definite colouring scheme.

As I have already pointed out, the colour modulation of the ribbons is given by the variable 'cm', at the beginning of 'ribbon.gnu', while the size of the figure with respect to the range of the values in 'ribbon.dat' is given by 'mult'. If you set it to something smaller, than 1, be prepared for some clipping of the data!

Finally, these two scripts will produce the graphs automatically for you, the only thing you have to set is the name of the file, and, of course, you have to have the data in the format given at the very beginning of my previous post, i.e., the data points must be arranged in a matrix.

Tomorrow (or whenever I have time), I will show what else this script can be used for, with some small modification. So long!

Making a ribbon chart in gnuplot

2009-08-06T10:58:00.001-07:00

While I am not a vary big fan of this kind of chart, I admit that on occasion, it might lend a refreshing new look to a graph. Besides, we can use it for something else, but more on that later. So, we will set to make this graph

from a data file similar to this

0.857354 0.943516 0.546751 -0.588053 -0.548767
0.551386 0.896607 0.163013 -0.166411 -0.442178
0.243002 0.310189 0.465029 0.0790944 0.378762
-0.24 -0.0625021 0.226812 0.48102 0.451281
-0.422674 -0.444781 0.314025 0.626378 0.822247
-0.513352 -0.893642 0.220384 0.800136 1.1908
-0.660823 -0.513169 0.316238 1.08866 1.25814
-0.0962106 -0.130051 0.0965893 0.944557 1.01623
-0.0152513 -0.0945605 0.468029 0.240179 0.556154
0.489453 0.726084 0.398192 -0.162081 -0.167408
0.876023 0.997108 0.139124 -0.281058 -0.581462
0.945061 1.19238 0.359603 -0.706723 -0.687105
0.860648 0.907059 0.217838 -0.540168 -0.592025
0.901743 1.02898 0.242619 -0.337589 -0.584658

which we will call ribbon.dat. There is an easy way to make this graph, and this solution requires an external script. Tomorrow, I will show a solution, if a bit convoluted, entirely in gnuplot, so that even those of you, who do not fancy a gawk one-liner, can make a similar figure.

The basic idea is to take the columns of ribbon.dat one by one, and turn them into a surface plot, whose derivative in the x direction is zero, while it takes the successive values of the column along the y axis. The gawk one-liner that does this transformation for us is given here

#!/bin/bash
gawk -v c=$2 '{printf "0 %d %g\n1 %d %g\n\n", NR-1, $c, NR-1, $c}' $1

I.e, we take the column designated by the value of 'c', and instead of listing the members, we print out six times the length of the column numbers, of the form

0 0 0.857354
1 0 0.857354

0 1 0.551386
1 1 0.551386

0 2 0.243002
1 2 0.243002

0 3 -0.24
1 3 -0.24

0 4 -0.422674
1 4 -0.422674
...

(If you look carefully, you will notice that this is nothing but the first four elements in the first column.) In this list, the first column is the 'x' value, the second column is the 'y' value, and the third is the 'z' value. We repeat all 'z' values twice, once with 0, and once with 1 'x' value. By plotting this, we will produce a ribbon of width 1. All that remains is to call the script multiple times, and change the colouring accordingly. With these in mind, our gnuplot script reads as

reset
unset key; unset colorbox
set tics out nomirror
set border 1+2+4+8+16+32+64+256+512 back
a=0.0
b=0.3
C=5.0

set xrange [-0.5:C-0.5]; set yrange [0:14]; set zrange [-2:2]; set cbrange [-2:2]
C=C+1.0
set xlabel 'x axis [a.u.]'; set ylabel 'y axis [a.u.]';
set ticslevel 0

set table 'bg.dat'
splot x
unset table

set multiplot
set palette model RGB defined (0 0 0 0.8, 1 0.9 0.9 0.9)
splot 'bg.dat' u 1:2:(-2):3 w pm3d, '' u 1:(14):($2/14*4-2):3 w pm3d, '' u (-0.5):2:($1-2):(-0.5) w pm3d

unset border; unset tics; unset xlabel; unset ylabel; unset zlabel
set palette model HSV function 1-a/C, 1-gray, 1-a/C
sp '<./ribbon.sh ribbon.dat 1' u (a+b*$1):2:3 w pm3d
a=a+1
set palette model HSV function 1-a/C, 1-gray, 1-a/C
sp '<./ribbon.sh ribbon.dat 2' u (a+b*$1):2:3 w pm3d
a=a+1
set palette model HSV function 1-a/C, 1-gray, 1-a/C
sp '<./ribbon.sh ribbon.dat 3' u (a+b*$1):2:3 w pm3d
a=a+1
set palette model HSV function 1-a/C, 1-gray, 1-a/C
sp '<./ribbon.sh ribbon.dat 4' u (a+b*$1):2:3 w pm3d
a=a+1
set palette model HSV function 1-a/C, 1-gray, 1-a/C
sp '<./ribbon.sh ribbon.dat 5' u (a+b*$1):2:3 w pm3d
unset multiplot

Here, 'b' is the width of the ribbons, 'C' is the number of ribbons we want to plot, and 'a' is just a variable for convenience. If you do not want to have the fancy background, you can skip lines starting with 'set table' to 'unset border'. Of course, if you do this, then you have got to move the line 'unset border;...' to after the first plot. In the rest of the script, we plot the columns, one by one, and also change the colour palette, so that the ribbons will have different colours. You can skip those lines, too, if you are satisfied with one colour for all the ribbons. I should also point out that the repeated call to the script can be made automatic, if you put those three lines into another gnuplot script, and call that script via reread. I discussed how to this this in my last two posts, so if you are interested in that, you should also read those. In those posts, you will also find hints as to how to avoid having to specify the various ranges by hand, and how to get the required values with the help of gnuplot. It should not be hard to implement those things in this script, so it is easy to make a "self-running" ribbon-plot-maker.

Well, this was about the easy way to make ribbon charts. In my next post, I will show how to do this in gnuplot, without having to rely on an external script. Till then!

Plotting a matrix with bargraphs

2009-08-02T11:35:00.000-07:00

In some cases, it is common practice to plot a matrix with bargraphs, instead of a colour map, or something similar. For instance, in physics, people usually plot the density matrix (if it does not contain too many elements) in this way. I have already discussed how bargraphs can be generated, but that method called an external gawk script, which had the advantage that we hadn't got to know the number of elements in advance. Yesterday, I showed how we can avoid the use of an external script, provided we know how many points we want to plot. But could we make something better, and produce a plot entirely with gnuplot, and without knowing how many, and of what magnitude the elements are? We could, by making use of some gnuplot variables. This is what I am going to discuss today. Since I will build upon the tricks that we exploited yesterday, I would recommend that you read that post, if you haven't already done so.

Let us suppose that we have the following data file, called 'bar_matrix.dat', with 4 by 4 real numbers

1 0.5 2.2 0.5
1.2 0.6 -2.4 0.2
0.1 0.5 -1.8 -1.5
1.2 0.3 0.6 1.9

and we want to generate a plot similar to this:

Following the recipe from yesterday, we have to read out the numbers in the matrix, one by one, and plot the cuboids one by one. In order to give the impression of a 3D object, the three visible sides of the cuboids must be painted with different shades of the same colour, but that should not be a problem, we have already seen how to do that. However, there is one difficulty here: we will have to use multiplot, which means that in order to keep the relative size of the cuboids, we have to fix all three ranges, xrange, yrange, and zrange. In addition, determining the xrange and yrange serves not only some aesthetic purpose, but they will also control our 'for' cycle. So, what can we do to crack this nut?

In order to answer the question, we have got to think about how gnuplot plots: first, if no ranges are specified, data are read in, then based on certain values, the ranges are calculated. So, if we issue

plot 'something.dat' using 1:2

the first, and second columns of 'something.dat' will be read, and based on the minimum and maximum values of the first column, the xrange is calculated, and likewise for the yrange. But then these values are stored as a special variable in gnuplot. So, if you are interested in the minimum and maximum of your xrange, you can do

print GPVAL_DATA_Y_MIN
print GPVAL_DATA_Y_MAX

You can also print out all variables, both user-defined and internal, by typing

show variables all

Now, this means that we will have access to some properties of our data file, provided that the data file has been plotted at least once. This is what we will use in the following script, called 'bar_matrix.gnu'.

reset
filename="bar_matrix.dat"
w=0.4
FIT_LIMIT=1e-8
mult=1.05; cd=1.4

f(x,a)=(abs(x-a)<0.5?d:0)
R(x)=abs(2*x-0.5); G(x)=sin(x*pi); B(x)=cos(x*pi/2.0)

set view 60, 20
set parametric; set isosample 2, 2
unset border; unset tics; unset key; unset colorbox
set ticslevel 0
set urange [0:1]; set vrange [0:1]

splot 'bar_matrix.dat' mat
set zrange [mult*GPVAL_DATA_Z_MIN:mult*GPVAL_DATA_Z_MAX]
Y=GPVAL_DATA_Y_MAX+1
X=GPVAL_DATA_X_MAX+1

set xrange [1-w:X+2*w]

set multiplot

C=0; xx=1; yy=Y
call 'bar_matrix_r.gnu'

set palette model RGB functions 0.9, 0.9,0.95
splot -w+u*(X+3*w), -w+v*(Y+w), 0 w pm3d

C=1; xx=1; yy=Y
call 'bar_matrix_r.gnu'
unset multiplot

In addition to this file, we will need two more. Remember, in order to plot a matrix, we have to use two 'for' cycles, each requiring an external gnuplot script. So, we also have 'bar_matrix_r.gnu'

yy=Y
call 'bar_matrix_r2.gnu'
xx=xx+1
if(xx<X+1) reread

and 'bar_matrix_r2.gnu'

unset parametric
yy=yy-1
d=0.3
set yrange [*:*]
fit [0:Y] f(x,yy) filename u 0:xx via d
set yrange [-w:Y+2*w]
set parametric
r=R(yy/Y); g=G(yy/Y); b=B(yy/Y)
if(C==0 && d<0) \
set palette defined (0 r g b, 1 r g b);\
splot w*u+xx, w*v+yy, d w pm3d;\
r=r/cd; g=g/cd; b=b/cd;\
set palette defined (0 r g b, 1 r g b);\
splot w*u+xx, yy, v*d w pm3d;\
r=r/cd; g=g/cd; b=b/cd;\
set palette defined (0 r g b, 1 r g b);\
splot w+xx, yy+w*u, v*d w pm3d;

if(C==1 && d>0) \
set palette defined (0 r g b, 1 r g b);\
splot w*u+xx, yy+w*v, d w pm3d;\
r=r/cd; g=g/cd; b=b/cd;\
set palette defined (0 r g b, 1 r g b);\
splot w*u+xx, yy, v*d w pm3d;\
r=r/cd; g=g/cd; b=b/cd;\
set palette defined (0 r g b, 1 r g b);\
splot w+xx, yy+w*u, v*d w pm3d;

if(yy>0) reread

Having seen the scripts, let us discuss what the three files do. The beginning of 'bar_matrix.gnu' should be familiar by now, there is nothing new in it. The first interesting thing happens where we do a dummy plot of the matrix, just to extract the number of rows and columns, and the range of the data values. We use GPVAL_DATA_Z_MIN and GPVAL_DATA_Z_MAX to set the common zrange of all consequent figures, and GPVAL_DATA_X_MAX and GPVAL_DATA_Y_MAX to determine the number of iterations (and the xrange, and yrange, of course). Note that we use a multiplier, 'mult', defined at the beginning of the file, so that the top and bottom of the plot will definitely not be clipped. Then we clear our plot by setting the multiplot, and calling 'bar_matrix_r.gnu'. This file controls the 'for' cycles over rows. It does nothing, but re-sets the value of the inner 'for' cycle, increments its own variable, and calls 'bar_matrix_r2.gnu'. In this second 'for' loop is where the actual plotting takes place. In that file, we call our fitting routine, set the palette, and plot the 3 sides of the cuboids. Note that we want to paint the sides in different colours, so after plotting the top, we produce a darker shade by dividing all RGB values by 'cd', and after plotting the front side, we reduce the RGB values once more, to draw the right hand side.

Now, this script contains an 'if' statement, with the variable 'C' and 'd'. The reason for this is that we call 'bar_matrix_r.gnu' twice: first we plot all bars that represent a negative value, then draw the z=0 plane, and finally plot the positive values. This is what happens in the last couple of lines of 'bar_matrix.gnu'. If you are certain that your matrix contains positive values only, you can get rid of the 'if' statement, and also of the first call of 'bar_matrix_r.gnu'.
If you now call 'bar_matrix.gnu' from gnuplot, you will get the figure at the beginning of this plot. You can then add labels and so on, as you wish. I would also like to point out that human intervention is required only at the very beginning of 'bar_matrix.gnu', where we specify the name of the file that we want to plot. Otherwise, everything is automatic, and will just happen by the stroke of the key. You should also keep in mind that we had a dummy plot, so if you want to plot into a file, you have got to set the terminal after this plot, lest you should not have the dummy plot in your file.

Pie charts - entirely in gnuplot!!!

2009-08-01T14:48:00.000-07:00

Some time ago, I discussed various ways to produce a pie chart in gnuplot. I have looked at the search terms that bring people to my blog, and it seems that this topic is one of the more popular, so I thought that it might be worthwhile to explore the issue more.

Since I am a Linux fan, I am quite fine with calling various scripts from gnuplot: pipes are rather convenient, and I can write a small script which does the data processing. But I see the downside of it, too: these solutions (while you can make it work under windows) will require extra steps, if not run under linux. I also understand that you might not want to delve into the nuances of gawk, for instance. So, I was wondering whether we could do everything in gnuplot, without relying on something external. And the answer is, yes, we can! (Otherwise, I wouldn't be writing this post at all...) What we are going to do is probably the dirtiest of hacks, for we will use a function of gnuplot, which was never intended to be used in such a ramshackle way. But we want to produce graphs, and this is not a coding beauty contest, after all!

If you recall, we had a data file with two columns: one designating years, the other one containing the values belonging to those particular years. Something like this, which we will call pie.dat:

You might also recall that the way in which we produced the pie chart was to plot arcs of a circle, the parameters determined by the second column. Now, the problem was (and this is why we had to use an foreign script) that those parameters cannot be set at run-time, so to speak, we had to hard-wire them into the gnuplot script. So, the question is really how we can access individual values of a data file, say, the 5th number in the 2nd column, and then do this repeatedly. The snag is that gnuplot hasn't a dedicated function to perform this task, so we have to look for a function that is dedicated to something else. Since there are not too many gnuplot functions that operate on a file, we can easily find the one that "returns" a value. We will use the fit function, and use the fitting parameter as a means of returning the sought-after value from the file. I know, I know! This is the ultimate abuse of 'fit', and we shouldn't do this at all! But what the heck! You haven't got to tell anyone how you produced that bloody figure!

Now, we have to find a proper fitting function. Remember, we want to pick the value of a particular number in the second column, say. Well, being a physicist, I would say that the Dirac-delta will do the job. However, we have to soften our stance a bit, and use something that is gnuplot-friendlier. For better or worse, I will choose the following function

f(x,a)=(x>a-0.5?(x<a+0.5?b:0):0)

You can recognise our old friend, the ternary operator, making its appearance twice in the definition of f(x,a). So, we first check, if the value 'x' is larger, than a-0.5. If so, we check whether it is smaller, than a+0.5. If this condition is fulfilled, we assign the value 'b', otherwise 0. By plotting it, you can see it for yourself that this function is a rectangle of height 'b' and width 1, centred on 'a', and 0 everywhere else. 'b' is going to be our fit parameter. Now, if you fit this function as

fit [0:7] f(x,2) 'pie.dat' u 0:2 via b
print b

then 0.2 will be printed. But that is the value of the 3rd number in the second column of pie.dat! (The rows are numbered starting with 0, that is why 2 means the 3rd row.) Since the value 'b' can be used in any subsequent expressions, definitions etc. as a number, we have found a way to extract any single number from any data file.

In order to proceed, we have to find a method to step through the rows of a column, one by one. For this purpose, we will use the reread command of gnuplot. You can learn the basic idea by issuing ?reread and ?if. 'reread' just repeatedly reads the file invoked by the last load command, and we can use 'if' to set some criterion as to how many times this repeated loading should take place. This is a primitive 'for' cycle, but it will do. (In gnuplot 4.3, the option 'for' was introduced in the 'plot' command, but that would not do too much good here, for we still have got to figure out the actual numbers.)

With these in mind, we write the following two scripts, the first of which named pie.gnu

reset
os=1.3
FIT_LIMIT=1e-8
L=6.0
f(x,a)=(x>a-0.5?(x<a+0.5?b:0):0)
r(x)=abs(2*x-0.5); g(x)=sin(x*pi); b(x)=cos(x*pi/2.0)

set view 30, 20
set parametric
set isosample 2, 2
unset border
unset tics
unset key
set ticslevel 0
unset colorbox
set urange [0:1]
set vrange [0:1]
set xrange [-2:2]
set yrange [-2:2]
set zrange [0:3]

A=0.0; D=0.0
set multiplot
# First, we draw the 'box' around the plotting volume
set palette model RGB functions 0.9, 0.9,0.95
splot -2+4*u, -2+4*v, 0 w pm3d
set palette model RGB functions 0.8, 0.8, 0.85

splot cos(u*2*pi)*v, sin(u*2*pi)*v, 0 w pm3d

call 'pie_r.gnu'
unset multiplot

while the second one named 'pie_r.gnu'

unset parametric
b=0.3
set yrange [*:*]
fit [0:L] f(x,D) 'pie.dat' u 0:2 via b
B=b
fit [0:L] f(x,D) 'pie.dat' u 0:1 via b
D=D+1.0
set palette model RGB functions r(D/L), g(D/L), b(D/L)
set parametric
set yrange [-2:2]
set urange [A:A+B]
set label 1 "%g", b at os*cos(2*pi*(A+B/2.0)), os*sin(2*pi*(A+B/2.0)), 0.2 cent
splot cos(u*2*pi)*v, sin(u*2*pi)*v, 0.2 w pm3d
A=A+B
if(D<L) reread

Now, let us see what is happening here. In the first script, 'os' will be used to place the labels later on. The gnuplot variable FIT_LIMIT is set to that value, so that the fit values are more accurate. It might be necessary to change it, if your labels are not what you expect. For the arcs, it should not really matter, because the default value is going to be accurate enough for any practical plots. 'L' is the number of data that we have (note that data are numbered starting with 0). Next, we define 4 functions. The first one is our fit function, the other three are used to colour the arcs. You can change these, if you are not satisfied with the colour scheme that you get. Any three functions will do, which are defined on [0:1], and return with a value in [0:1]. You can read about this in the post in which I discussed phonged surfaces, sometime in late May.

The next couple of lines set up the various ranges of our plot. I wrote about this in my first pie chart post, and the lines giving the background of the plot should also be familiar from that post. 'A' and 'D' are our control variables that we manipulate in 'pie_r.gnu'. We then draw the shadow of the pie, and finally, call 'pie_r.gnu'.

Let us take a look at 'pie_r.gnu'. First, we extract the value in the second column, and then in the first one. We will use this latter one to produce the label. Note that we have got to unset the parametric plot, otherwise, the fitting function would not work. Also note that we re-set the yrange. This is necessary, because the actual plot is in the [-2:2] range, while the fit is on values around 1990. Then we increment the value 'D' (this is the ordinal number of the row that we are currently processing), and re-set our palette, using the three functions, r(x), g(x), and b(x) that we defined in 'pie.gnu'. Then, using the extracted values, we set the range of the parametric plot, and define the label, and plot the arc. Finally, we check, if we have called the script enough times. Loading 'pie.gnu' will produce the following graph:

A couple of comments about this script: as I already mentioned, you can change the colour scheme by using various functions. However, if you are really lazy, you can simply generate three random numbers by

r=rand(0); g=rand(0); b=rand(0)

and assign these values to the next palette.

Second, you can easily implement an "explosion", by shifting one or several of the arcs by shift*cos(2*pi*(A+B/2)) and shift*sin(2*pi*(A+B/2)), based on some condition that you set. It should not be hard, either, to plot a real 3D pie char based on the scheme that I outlined about two months ago.

Third, if the sum of your numbers is not normalised to one, you can easily fix it by adding an extra loop to your script: if you have a script called 'pie_r2.dat', containing

unset parametric
b=0.3
set yrange [*:*]
fit [0:L] f(x,D) 'pie.dat' u 0:2 via b
D=D+1.0
G=G+b
if(D<L) reread

and call this in 'pie.gnu' immediately before 'pie_r.gnu', then 'G' will just be the sum of all numbers in column 2, and you can use this to normalise the numbers when you call 'pie_r.gnu'.

And last, the only thing you need in advance is the number of records you want to process, 'L'. This is the only thing you have to set by hand, all the rest is automatic. If you want to learn how to avoid this small "difficulty", you should read the post on the 2nd of August.

I hope that this script can convince people that we haven't got to rely on any outside tool, and can still produce a decent pie chart in gnuplot.I also feel that this approach is, in some sense, much cleaner, than the one that was based on gawk. Some people will probably contest this statement, but doing everything in gnuplot ensures that we retain the platform independence.

Adding shadow to the key

2009-07-29T14:02:00.000-07:00

I think, it changes the impression a graph can make by quite a lot, if it looks sort of three-dimensional. Even a tiny twist to the boring 2D lookout can make a difference, though it doesn't give any new information to the reader. However, in a presentation, it is almost expected that one produces "exciting" graphs. Today, we will discuss an easy way to add a slight shadow to the key on a graph, as if it were lifted from the plane of the curves. At the end of our exercise, we will have the following figure

I believe the steps are quite straightforward, so I won't spend too much time on explaining every detail. Here is the script that we need

reset

xl=0; xh=1; yl=-1; yh=1;
eps=0.01;
rx=0.6; ry=0.8; kw=0.35; kh=0.15
lh=0.06; al=0.1

key1="First function"
key2="Second function"

set table 'shadowkey.dat'
splot [xl:xh] [yl:yh] x/(xh-xl)
unset table

set object 1 rect from graph rx,ry rto kw,kh fc rgb "#aaaaaa" fs solid 1.0 front lw 0
set object 2 rect from graph rx-eps,ry+eps rto kw,kh front fs empty
set label 1 at graph 1.1*al+rx, ry+2*lh key1 front
set label 2 at graph 1.1*al+rx, ry+lh key2 front
set arrow from graph rx, ry+2*lh rto al, 0 lt 1 lw 1.5 nohead front
set arrow from graph rx, ry+lh rto al, 0 lt 3 lw 1.5 nohead front

unset colorbox
unset key
set palette defined (0 "#8888ff", 1 "#ffffff")
plot [xl:xh] [yl:yh] 'shadowkey.dat' w ima, \
x*x*exp(-x) lw 1.5, cos(13*x)*exp(-3*x) lt 3 lw 1.5

I collected all variables at the beginning of the script, so that it will be easier to adapt it to any situations. The first four numbers define the xrange and the yrange. The next four numbers specify the position and the size of the key. We will have to make our own key, and the size of the key will depend on the particular terminal one uses. This is why it is handy to define them at the beginning of the script, so if you are not satisfied with the results, you can easily adjust both the position and the size. The next number (lh) gives the hight of one key line, while 'al' will be the length of the line which represents the curve. 'key1' and 'key2' are just two arbitrary strings, holding the text of the key.

The next three lines are necessary only if you want to draw the background gradient, but if you are happy with a white plot, you can skip this. I should also point out that the first four numbers that are needed for the xrange and yrange, are needed only because we have to "synchronise" the plot with its background. If you don't want the gradient, you can drop these definitions, and you haven't got to specify the range in the plot either. This might be useful, when you don't know the plotting range beforehand, and want to let gnuplot calculate it for you.

The next step is to draw two rectangles in the front of the graph: one gray, and one white. Note that the rectangle in gray has no boundary, i.e., that is set to zero width. Also note that we specify the coordinates in terms of the numbers we defined at the beginning, so both rectangles will change accordingly, if you change those numbers. This means that the white rectangle and its shadow will always be linked, given by the variable 'eps'. When we are done with the white rectangle, we can place the text of the keys, and the two lines representing our plots. For this purpose we draw two arrows without heads, and with the linetypes that we will use for plotting the curves.

The final step is to actually draw the curves. Before that, we plot our file, 'shadowkey.dat', so the graph will have a background gradient. If you skipped those three lines writing the data file to disc, you should replace the plotting command by

plot [xl:xh] [yl:yh] x*x*exp(-x) lw 1.5, cos(13*x)*exp(-3*x) lt 3 lw 1.5

It should go without saying that in this case, you can also drop the palette definition, for it will not be used. Well, this is it for today. I know that this was a simple trick, but we can't have something complicated every day!

Maps - Contour plots with labels

2009-07-26T13:01:00.000-07:00

So, you have always wondered how on Earth one can make a real map with gnuplot. Well, there is a simple and a not-so-simple way to this. First, let us see the simple way. Since it is simple, this method won't have one of the features, isoline labelling, that the second one has. As we go along, the map will become more and more complicated, but I hope that I set the right pace, and it will be easy to follow.

A map is nothing but a colour-coded 3D plot, with the isolines attached to it. We will produce a simple map using the function

sin(1.3*x)*cos(0.9*y)+cos(.8*x)*sin(1.9*y)+cos(y*.2*x)

I don't think that this function has any particular meaning, but it looked quite all right, at least, as far as maps and isolines are concerned. If you have a matrix to plot, you can replace this function with that file. In principle, we could plot both the colour-coded map and the isolines from gnuplot, but we will have much greater flexibility, if we first direct the plots to a file, and then call the data from those files. One of the advantages of doing this is that in this way, we can make sure that the various plots have the same size. This requires a small overhead in terms of scripting, namely, we have got to issue the command

set table 'something.dat'
splot something
unset table

While this might appear superfluous, there are good reasons to do this. If we plot our function through the file 'something.dat', then we can use the 'with image' modifier of the plot command, and this means that the plot will be of the same size as would a normal 2D plot. Otherwise, if we use

set pm3d map
splot something with pm3d

the plot will actually be a bit smaller. The reason behind this is that in a 3D plot, we have to have some space for the 'z' axis, and even if we drop it in the map view, the space-holder for the 'z' axis is still there, therefore, gnuplot makes the whole plot a bit smaller. If, however, we plot the map through a file, we can use 'plot', in splot's stead, therefore, the 'z' axis will not appear anywhere in the processing of the plot.

After this interlude, let us see our first version of the map.

reset
f(x,y)=sin(1.3*x)*cos(.9*y)+cos(.8*x)*sin(1.9*y)+cos(y*.2*x)
set xrange [-5:5]
set yrange [-5:5]
set isosample 250, 250
set table 'test.dat'
splot f(x,y)
unset table

set contour base
set cntrparam level incremental -3, 0.5, 3
unset surface
set table 'cont.dat'
splot f(x,y)
unset table

reset
set xrange [-5:5]
set yrange [-5:5]
unset key
set palette rgbformulae 33,13,10
p 'test.dat' with image, 'cont.dat' w l lt -1 lw 1.5

In the first couple of lines, we plot the actual function into a file, called 'test.dat'. There will be some 25000 data points in this file, for there are 100 samples (default), and 250 isosamples. Having written these points into a file, we plot the contour. We just happen to know that the value of f(x,y) is between -3 and 3, so we set the contour levels at -3, -2.5...2.5, 3. When we are done with all this, we simply reset our plot, set the x and yrange, specify the palette that we want to use with the colouring, and call the plot command. For the actual data points, we use the 'with image' modifier, while for the contours, we increase the linewidth to 1.5, instead of using the default value of 1. We, thus, have just produced the following image.

Now, this was sort of standard, but the question inevitably comes up, whether we could do something more with this. When one has contour lines, one wants to label them, I presume. Petr Mikulik has a gawk script that will do this, in particular, it will place the appropriate labels next to the contour lines. I will discuss another approach here. The main difference with Petr Mikulik's method is that here we will put the labels on the contour lines, with some white space, of course.

First, it is quite instructive to look at the file containing the contour lines. In order not to cobble the plot, we will restrict it to the range [-5:0] (x) and [2:5], but it is really just for the sake of simplicity. With this modification, the file should read as

#Surface 0 of 1 surfaces

# Contour 0, label:        2
-0.482969  3.59036  2
-0.505051  3.58024  2
-0.509852  3.57831  2
-0.539895  3.56627  2
-0.555556  3.55994  2
-0.572136  3.55422  2
...

# Contour 1, label:      1.5
-1.11111  4.57693  1.5
-1.10874  4.57831  1.5
-1.0877  4.59036  1.5
-1.06617  4.60241  1.5
-1.06061  4.60546  1.5
-1.04273  4.61446  1.5
...

What this means is that the 0th contour line is drawn at z=2, and then the (x,y,z) values are listed. Of course, since we know that this particular line is at z=2, the last coordinate doesn't play any role, but it is listed all the same. When we plot this file, consecutive (x,y) points will be connected by a straight line segment. The next contour is at 1.5 etc. It might happen that one contour line is made up of several blocks, both for technical, and fundamental reasons. The technical reason is gnuplot's inner procedure of computing the an isoline, while the fundamental is simply that there might be several disconnected regions with the same isolevel. Anyway, if there are several blocks, they are separated by a blank line within the same contour.

In order to create the labels and the white spaces for them, we will use the following short script:

#!/bin/bash

gawk -v d=$2 -v w=$3 -v os=$4 'function abs(x) { return (x>=0?x:-x) }
    {
            if($0~/# Contour/) nr=0
            if(nr==int(os+w/2) && d==0) {i++; a[i]=$1; b[i]=$2; c[i]=$3;}
            if(abs(nr-os-w/2)>w/2 && d==1) print $0
            nr++
    }
    END {   if(d==0) {
                    for(j=1;j<=i;j++)
                    printf "set label %d \"%g\" at %g, %g centre front\n", j, c[j], a[j], b[j]
            }
    }' $1

This script operates on the data file of contour lines that we had before, and simply takes out a couple of points of the contours, while keeping one of the coordinates as the position of the labels. You can change the length of the white space by setting the 'w' argument, while the position of the white space can be modified by changing the 'os' argument. The first argument, 'd', determines what we want to do with the script: make the labels, or tweak the input file.
Now, we can call the script from our gnuplot script as follows:

reset
set xrange [-5:0]
set yrange [2:5]
set isosample 150, 150
set table 'test.dat'
splot sin(1.3*x)*cos(.9*y)+cos(.8*x)*sin(1.9*y)+cos(y*.2*x)
unset table
set cont base
set cntrparam level incremental -3, 0.5, 3
unset surf
set table 'cont.dat'
splot sin(1.3*x)*cos(0.9*y)+cos(.8*x)*sin(1.9*y)+cos(y*.2*x)
unset table
reset
set xrange [-5:0]
set yrange [2:5]
unset key
set palette rgbformulae 33,13,10
l '<./cont.sh cont.dat 0 15 0'
p 'test.dat' w ima, '<./cont.sh cont.dat 1 15 0' w l lt -1 lw 1.5

There are two more twists that we could add to the figures above. One is that the labels can be rotated in such a way that they are parallel to the curves. This we can easily achieve by adding two more lines to our script above: we take the first and last dropped coordinate, and calculate the angle that that line segment would make with the horizontal. Using this angle, we rotate the labels accordingly. With this modification, the script now becomes

#!/bin/bash

gawk -v d=$2 -v w=$3 -v os=$4 'function abs(x) { return (x>=0?x:-x) }
    {
            if($0~/# Contour/) nr=0
            if(nr==int(os+w/2) && d==0) {a[i]=$1; b[i]=$2; c[i]=$3;}
            if(nr==int(os+w/2)-1 && d==0) {i++; x = $1; y = $2;}
            if(nr==int(os+w/2)+1 && d==0) r[i]= 180.0*atan2(y-$2, x-$1)/3.14
            if(abs(nr-os-w/2)>w/2 && d==1) print $0
            nr++
    }
    END {   if(d==0) {
                    for(j=1;j<=i;j++)
                    printf "set label %d \"%g\" at %g, %g centre front rotate by %d\n", j, c[j], a[j], b[j], r[j]
            }
    }' $1

and the resulting graph is shown here:

Note that the only difference between this and the previous script is the following two lines


if(nr==int(os+w/2)-1 && d==0) {i++; x = $1; y = $2;}
if(nr==int(os+w/2)+1 && d==0) r[i]= 180.0*atan2(y-$2, x-$1)/3.14

where we calculate the appropriate angles.

The second twist is that it might make interpretation of the figure easier, if the contours are coloured, as are the corresponding labels. We can use the script above with some small modification. The key is to use the 'index' modifier. Recall that the contours are rendered in blocks, and with the help of the 'index' keyword, we can specify which block we want to plot. We, therefore, modify our script as follows:

#!/bin/bash

gawk -v d=$2 -v w=$3 -v os=$4 'function abs(x) { return (x>=0?x:-x) }
   {
           if($0~/# Contour/) nr=0
           if(nr==int(os+w/2) && (d%2)==0) {a[i]=$1; b[i]=$2; c[i]=$3;}
           if(nr==int(os+w/2)-1 && (d%2)==0) {i++; x = $1; y = $2;}
           if(nr==int(os+w/2)+1 && (d%2)==0) r[i]= 180.0*atan2(y-$2, x-$1)
           if(abs(nr-os-w/2)>w/2 && (d%2)==1) print $0
           nr++
   }
   END {   if(d==0) {
                   for(j=1;j<=i;j++)
                   printf "set label %d \"%g\" at %g, %g centre front rotate by %d tc lt %d \n", j, c[j], a[j], b[j], r[j], j

           }
           if(d==2) {
                   printf "plot \"test.dat\" w ima, \"cont.plt\" index 0 w l lt 1,\\\n"
                   for(j=2;j<i;j++) printf "\"\" index %d w l lt %d,\\\n", j-1, j
                   printf "\"\" index %d w l lt %d\n", i-1, i

           }
   }' $1

Calling the script produces the following graph:

Note that we had to call plot through a new temporary file, cont.plt, which should be produced by redirecting the output of

cont5.sh 2 15 0 > cont.plt

Grid with multiple colours

2009-07-15T13:32:00.000-07:00

A couple of days ago, I very badly needed a graph on which all four axes were used, and in addition to that, I needed a grid on the graph. Now, the problem is, if you just plot a grid, it will be hard to follow which grid line corresponds to which axis. So, I decided to colour both the grid, and the axes, and use the same colour for corresponding grid lines and axes. Well, this is what I decided, but then it turned out to be a rather hard nut to crack. (In gnuplot 4.3, the situation is a bit better, I will discuss that tomorrow.) I have found a work-around, however, so if you are interested, you can find the details below. But let us see the graph first!

The standard way of defining a grid is a line similar to this:

set grid ytics lt 0 lw 1 lc rgb "#880000"

First, we define where we want to have the grid (everywhere, where there is an ytics), then the linetype, the line width, perhaps, and the colour. The problem is that when you issue this command multiple times, for some funny reason, only the last colour setting will take effect, i.e., if after the line above, you happened to have

set grid y2tics lt 0 lw 1 lc rgb "#008800"

then gnuplot would place the grid lines at y2tics all right, but the grid lines at ytics would also be green, and this is not what we want. The way out of this difficulty is to trick gnuplot into thinking that we have more than one plot, and re-set the grid before each new plot. Of course, we don't actually want to plot anything after the first one, therefore, we will just plot 1/0. It is a vicious way to deceive gnuplot, but it does the trick. Now, here is the script

reset
h=4.136e-15 # eV s
hbar=6.57e-16 # eV s
c=3e17  # nm/s

set multiplot
unset key
set xlabel 'Wavelength [nm]' tc rgb "#0000ff"
set x2label 'Temperature [K]'
set ylabel 'Energy [eV]' tc rgb "#880000"
set y2label 'dE/d{/Symbol l} [meV/nm]' tc rgb "#008800"
set xrange [200:2000]
set x2range [0:10]
set yrange [0:6]
set y2range [1:-7]
set y2tics 1,-1,-7
set ytics nomirror tc rgb "#880000"
set y2tics tc rgb "#008800"
set xtics nomirror tc rgb "#0000ff"
set x2tics
set mxtics 2
set mytics 2
set grid xtics lt 0 lw 1 lc rgb "#0000ff"

plot h*c/x w l lc rgb "#880000", -1000*h*c/x/x axes x1y2 w l lc rgb "#008800"

unset grid
set xlabel " "
set xtics format " "
set x2tics format " "
set x2label " "
set ylabel " "
set y2label " "
set grid x2tics lt 0 lw 1 lc rgb "#000000"
plot 1/0

unset grid
set grid ytics lt 0 lw 1 lc rgb "#880000"
plot 1/0

unset grid
set grid y2tics lt 0 lw 1 lc rgb "#008800"
plot 1/0

unset multiplot

First, we define a couple of constants, and then set up the axis labels. Note that when doing so, we assign a colour to each axis, so that we will know which axis the plot belongs to. Of course, we could do this with arrows, but since we want to have coloured grids, it is better to do it in this way. Next, we set the various ranges, and then set up the first grid line, which will be drawn at every xtics, and in blue. (Blue was the label colour of x.) After plotting the curves, we unset the grid, and we also take off the labels and ticlabels. Note that when doing so, we, in fact, re-set the labels, but with an empty string. The reason for this is that, if we simply take off the labels, then the margins of the graph will change, and when we plot our next graph (which is empty, but we still need the grid), it would be of a different size, and would not align with our original plot. Having re-labelled the graph, we re-set the grid, this time putting the lines on x2tics in black, then on ytics in a dark shade of red, and finally, on y2tics, in dark green. At the very end, we unset the multiplot, i.e., escape from the whole plot.
This is a bit dirty, but if you are pressed for it, it should still be OK. Next time, I will show the easy route, provided you have gnuplot 4.3.

Bubble graphs with a different method

2009-06-28T04:45:00.000-07:00

Yesterday, I showed a simple way to producing bubble graphs with gnuplot. That method relied on manipulating the eps file. Today, we will try another option, this time plotting everything on the appropriate terminal. At the end, we will have the following figure:

First, let us have a look at the script, assuming that we have the following data to plot:

0 -0.0694726
0.20202 0.233484
0.40404 0.424311
0.606061 0.546688
0.808081 0.580043
1.0101 0.862214
1.21212 0.907601
1.41414 0.759692
1.61616 0.884879
1.81818 0.784566
2.0202 0.71774...

reset
f(x) = A*exp(-x*x/B/B)
rx=0.107071; ry=0.057876; A = 1; B = 0.2; C=0.5*rx; D=-0.4*ry
g(u,v) = (2*cos(u)*v*rx+C)*(2*cos(u)*v*rx+C)+(3.5*sin(u)*v*ry+D)*(3.5*sin(u)*v*ry+D)             
unset key; unset colorbox; set view map
set xrange [-0.15:5.2]; set yrange [-0.7:0.95]
set parametric; set urange [0:2*pi]; set vrange [0:1]                         
set isosamples 20, 20; set samples 30                                         
set palette model HSV functions 1, 1-f(gray), 1+2*f(gray)                     
splot cos(u)*rx*v+0.000000,sin(u)*ry*v+0.000000, g(u,v) w pm3d, \
cos(u)*rx*v+0.202020,sin(u)*ry*v+0.233484, g(u,v) w pm3d, \
cos(u)*rx*v+0.404040,sin(u)*ry*v+0.424311, g(u,v) w pm3d, \
cos(u)*rx*v+0.606061,sin(u)*ry*v+0.546688, g(u,v) w pm3d, \
cos(u)*rx*v+0.808081,sin(u)*ry*v+0.580043, g(u,v) w pm3d, \ ...

First, we define a Gaussian; this will be the colouring function, and then define a couple of variables that go into that function. Having done this, we define the function that determines the argument of f(x). The next couple of lines simply sets the ranges and the parametric plot with the parametric ranges, and finally, the number of samples. The last thing we have to define is the palette function. We choose red (i.e., the hue is equal to 1), and the saturation and value are given by 1-f(x) and 1+2*f(x), where the argument is the gray value. At this point, we are ready to plot the points in question. We will simply draw circles with origin x,y, where the x and y values are taken from the data file that we showed above. Note that in fact we are drawing ellipses with axes rx and ry. The reason for this is that the aspect ratio of the plot is not equal to 1, i.e., were we to draw circles, they would look ellipses on the plot. The value of rx and ry are determined by the plot ranges xrange and yrange. (We will see this in the gawk script below.) When plotting the points, we have to plot one circle for each data point, i.e., we have to call the plot function many times, while C and D give the centre of the white spot. increasing C or D will push the white points to the edge of the circles.

Now, a few words on the various parameters above. The value A determines how bright the bubble will be at its brightest point. 1 corresponds to white, values smaller than 1 give a darker tinge of red. B determines how tight the white spot is. Obviously, rx and ry are the size, so if you want to have smaller circles, you could scale them accordingly, keeping their ratio.

We could easily write a script that takes a data file with two (or more) columns, and turns it into a gnu script along the lines presented above. A possible implementation in gawk is here.

#!/bin/bash

gawk  '{
  if($0!~/#/) {
   x[i] = $1
   y[i] = $2
   if(i==0) { mx = x[i]; my = y[i] }
   if(i>0)   {
         if(max < x[i]) max = x[i]
         if(mix > x[i]) mix = x[i]
         if(may < y[i]) may = y[i]
         if(miy > y[i]) miy = y[i]
   }
   i++
 }
 }
 END { eps = 0.03
  lx = mix-eps*(max-mix)
  hx = max+eps*(max-mix)
  ly = miy-eps*(may-miy)
  hy =  may+eps*(may-miy)
  print "reset"
  print "f(x) = A*exp(-x*x/B/B)"
  printf "rx=%f; ry=%f; A = 1; B = 0.2; C=0.5*rx; D=-0.4*ry\n", 0.02*(hx-lx), 0.035*(hy-ly)
  print "g(u,v) = (2*cos(u)*v*rx+C)*(2*cos(u)*v*rx+C)+(3.5*sin(u)*v*ry+D)*(3.5*sin(u)*v*ry+D)"
  print "unset key; unset colorbox; set view map"
  printf "set xrange [%f:%f]; set yrange [%f:%f]\n", lx, hx, ly, hy
  print "set parametric; set urange [0:2*pi]; set vrange [0:1]"
  print "set isosamples 20, 20; set samples 30"
  print "set palette model HSV functions 1, 1-f(gray), 1+2*f(gray)"
  printf "splot "
    for(k=0;k<i-1;k++) {
         printf "cos(u)*rx*v+%f,sin(u)*ry*v+%f, g(u,v) w pm3d, \\\n", x[k], y[k]
    }
    printf "cos(u)*rx*v+%f,sin(u)*ry*v+%f, g(u,v) w pm3d\n", x[i-1], y[i-1]
 }' $1

At the beginning, we fill up the x[] and y[] vectors, while, at the same time, determining the minimum and maximum of these two vectors. Then we use these values to determine xrange and yrange (I defined them a little bit bigger than the minimum and maximum of the vectors, so that the circles are confined in the plot.), and the value of rx and ry, as well. At the end, we simply call the plot function with the arguments that we take from the x[] and y[] vectors.

I believe it should be fairly easy to modify the script, should you want to make some changes to it. Finally, a word of caution: since we use some 900 samples for each circle we plot, this is going to be reflected in the file size, if you use a vector output. As a comparison with the method that I discussed yesterday, the size of that postscript file was something around 20 kB, while the size of the file we would produce with the present method is about 600 kB. We have this big difference simply because yesterday we re-defined one of the symbols, while today we plot each point, without any reference to a particular symbol. Therefore, if you want to include the plot in a publication, it is better to use yesterday's method. For bitmap files, png, jpeg and the like, there should be no significant difference in size.

Making bubble graphs with gnuplot

2009-06-27T12:05:00.000-07:00

Today, I will try to explain how one can make a bubble graph with gnuplot. This is a low-level hack, and is probably not for the faint-hearted. On another occasion, I will discuss how this can be done in a more convenient way, at the price of producing larger files. (Well, this is true, only if one makes postscript figures, but for raster, like png, it shouldn't matter.)

One of gnuplot's most appealing features (at least, for me) is that the postscript file it produces is actually human-readable, and quite instructive. This makes it possible to easily change the plot even when it is in fact finished in the sense that gnuplot has produced it.
So, if you feel a bit adventurous, you could try to see what can be achieved by manipulating the postscript file. Plot the data points with point type 7 (as a matter of fact, it doesn't really matter which style you use, but in what follows, I will assume that we have point type 7), and redirect the output to a postscript file, i.e., you would have a gnu script as in

reset
unset key
set border back
set terminal postscript eps enhanced
set output 'bubble.eps'
plot 'bubble.dat' u 1:2 w p pt 7 ps 2

Simple enough, isn't it? (Do not forget to reset the terminal, once you are done with the plot. Otherwise, all your consequent graphs will be written to 'bubble.eps'.)

Now, that we have the postscript file, open it in any text editor, and search for the line


/CircleF { stroke [] 0 setdash hpt 0 360 arc }

This should be somewhere around line 161. This is just the definition of point type 7, and this definition is used in the actual plotting of your data, somewhere at the end of the postscript file. We will replace the line above by

/CircleF { stroke [] 0 setdash 0 0 0 arc
<< /ShadingType 3
/ColorSpace /DeviceRGB
/Coords [currentpoint exch hpt 1 mul sub exch hpt 0.2 mul add 0 currentpoint 150]
/Function << /FunctionType 2
/Function { 1 }
/Domain [0 1]
/C1 [ 0.8 0 0 ]
/C0 [ 1 1 1 ]
/N 1>> >> shfill} def

While I am not going to write a tutorial on the postscript language, we could certainly see what this does. What we have to keep in mind is that postscript works on stack variables, i.e., we can access various quantities in a top-to-bottom direction, with the last operand being at the top.

First, we take the current x and y coordinates (this is the first currentpoint instruction), and reverse the order of x and y (exch), subtract the value hpt from x (hpt 1 mul sub), reverse the order again, so we are back to the original x-y order, add 0.2*hpt to y (hpt 0.2 mul add). With this set of operations we defined the centre of the colour of our bubble, which is denoted by the next 0. (So, the centre of colour will be to the left and to the top of the geometric centre. North-west, that is.)

Then we define the origin as currentpoint, and the radius as 150. The next couple of lines define the function type and the domain of this function. /C1 and /C0 determine the colour of the object when the function in question takes on the value of 1 (this will be a bit dark shade of red) and 0 (which will be glowing white). The last line just fills the circle. Now, that we know how this works, we can tweak our figure, if needed. As I said, the size of the bubbles will be determined by the last number on the line beginning with /Coords, and we can move the centre of the colour by changing 0.2 to something else (this would be the y direction), or 1 (x direction). In this way, we can where the light appears to come from. Finally, the colour can be changed by replacing the lines in /C1 and /C0. All in all, we have produced the following figure. While this might appear a bit grainy, this is only an artefact of the postscript viewer. If you print it out, it's going to be all right.

As I promised at the beginning, I will show an easier route to achieving the same result. However, we can already see that since the manipulation of the postscript file doesn't require any prior knowledge about anything, this process can easily be scripted, so that you haven't got to replace the lines by hand every time you want to make bubbles.

Broken axis revisited

2009-06-26T13:43:00.000-07:00

In one of my first posts, I described a method by which one can draw broken axes in gnuplot. Now, that method suffers from many problems, most notably that it relied on gnuplot 4.3 (still in development), and one had to do almost everything by hand. Today we will look at a different method, which eliminates these problems, and which requires very little human intervention. So, we would like to have something like this graph:

(This is just a sin function, broken at 4.5, and for the right hand side, displaced by 3). The recipe that we are going to follow takes only a couple of lines. But, let us see the code!

reset
A=4.5
B=3.0
C=0
D=10
E=1
eps=0.05*E
eps2=0.005*(D-C)
unset key
f(x) = (x<A?1:1/0)
g(x) = (x>A?1:1/0)
h(x) = (x<A?x:x+B)
set xtics 0, 2, A
set xtics add (gprintf("%.0f", 6+B) 6)
set xtics add (gprintf("%.0f", 8+B) 8)
set xtics add (gprintf("%.0f", 10+B) 10)

set xlabel 'Time [s]'
set ylabel 'Position [m]'
set yrange [-E:E]
set arrow 1 from A-eps2, -E to A+eps2, -E nohead lc rgb "#ffffff" front
set arrow 2 from A-eps2, E to A+eps2, E nohead lc rgb "#ffffff" front
set arrow 3 from A-eps-eps2, -E-eps to A+eps-eps2, -E+eps nohead front
set arrow 4 from A-eps+eps2, -E-eps to A+eps+eps2, -E+eps nohead front
set arrow 5 from A-eps-eps2, E-eps to A+eps-eps2, E+eps nohead front
set arrow 6 from A-eps+eps2, E-eps to A+eps+eps2, E+eps nohead front
plot [C:D] f(x)*sin(x) w l lt 1, g(x)*sin(h(x)) w l lt 1

The first several lines are just definitions: 'A' is the position at which we want to break the line, 'B' is the value by which the right hand side of the graph will be shifted, 'C', 'D' and 'E' are just definitions for the x and y ranges, while 'eps' and 'eps2' will be needed for the drawing of the slanted tics representing the discontinuity in the axes.

The first really important definition is that of f(x), g(x) and h(x), which are helper functions (green lines). In all three cases, we use the ternary operator that I discussed in my previous post. Basically, f(x) returns 1, if x is smaller than 'A', g(x) does just the opposite, while h(x) produces a shift of its argument, if x is larger than 'A'.

Then we set the labels on the x axis (blue line in the code). But watch out: we do it only up to 'A', because the axis is going to be broken at 'A'. In order to make up for the missing tic marks, we add them in the next three lines (red). We can save some work, and can eliminate the possibility of messing something up, if we ask gnuplot to compute the values for us; this is done by the string conversion function gprintf, and adding the converted string.

set xtics add (gprintf("%.0f", 10+B) 10)

The last bit is drawing the break points on the axis. In order to do so, we draw 6 arrows with no heads. Note that the first two arrows are white, and they serve as a beauty plaster: we put them on top of the axes, to give the impression that the axes are broken. Also note that since arrows are drawn earlier than axes and borders, we have got to explicitly instruct gnuplot to bring the arrows to the front of the figure. The last four line segments are slightly slanted, to give a more appealing look to the graph. Once we are done with the set up of the figure, we can plot the function. The only thing we have to keep in mind is that the right hand side should be shifted by 'B', which is done by calling h(x) in the argument of sin(x).

The procedure for the y axis is very similar, so I just dodge the discussion of that. I should point out that you can easily change the angle of the last 4 arrows by multiplying the shift of the y coordinates of the end points by something, and it is just as easy to move them farther from each other, if you replace the definition of 'eps2'. Modifying 'eps' will change their length, but not their angle.

Till next time!

Plotting an inequality in 2D

2009-06-20T12:38:00.000-07:00

In one of the comments, someone asked how an inequality can be plotted in gnuplot. Well, there are two solutions (at least) for this problem. One of them is discussed in the following page, if you need only the curve separating the region of the x-y plane that fulfils the inequality from that that doesn't. But the question was how to shade these regions. We will look at the inequality

x³ - 2xy + y³>0

The way to do this is to define a function that has a constant value, if the inequality holds true, and has an indefinite value, if this is not the case. The value of 1/0 is quietly ignored in gnuplot, thus we can define our function as follows

f(x,y)=(x*x*x - 2*x*y + y*y*y>0)?1:1/0

where we used the ternary operator, which always has the form

A?B:C

i.e., returns 'B', if 'A' is true, and 'C' otherwise. (Incidentally, the value 'B' can be a ternary operator in itself, and in this way, multiple conditions can be imposed.) In the function above, f(x,y) = 1, if x³ - 2xy + y³>0, and f(x,y) = 1/0 (undefined), if x³ - 2xy + y³<0. Once we defined our function, we can plot it in the usual way. The complete script and the figure is given below.

reset
f(x,y)=(x*x*x - 2*x*y + y*y*y>0)?1:1/0
unset colorbox
set isosample 300, 300
set xlabel 'x'
set ylabel 'y'
set sample 300
set pm3d map
splot [-2:2] [-2:2] f(x,y) t 'x^3 - 2xy + y^3>0'

Stacking shiny histograms

2009-06-20T11:37:00.000-07:00

As I promised last week, we will re-visit the problem of the shiny histograms, and we will stack them. After that discussion, it will be plain sailing, because we have everything at our disposal, we have simply got to re-arrange the columns a bit. Obviously, there are two ways of stacking histograms, and I will give a script for both cases. Thus, we will achieve something like either of these plots. (The data file is the same as that in the post on 14th June.)

The difference of these two plots is just that the second one is normalised, so all values are divided by the sum in a particular group.
Now, as for the script generating the first graph, that reads as this:

reset
f(x, a, b, c) = a*exp(-(x-b)*(x-b)/c/c)                          
A = 0.6; B = 0.5; C = 0.2; D=0.5                                 
unset key; unset colorbox; unset xtics; set pm3d map             
set yrange [0:6.6]                                          
set parametric; set urange [0:1]; set vrange [0:1]               
set xrange [0:5.5]                                          
set multiplot                                                    
set isosamples 100, 100                                          
set palette model HSV functions 2.0/3.0, 0.4-0.4*gray, (4.0+2.0*gray)/6.0
splot 5.5*u, 6.6*v, 2*u+v w pm3d                       
set isosamples 2, 2                                              
set palette model HSV function 1.0, 1.0-f(gray, A, B, C), (1.0+2.0*f(gray, A, B, C))/3.0                                                         
splot D*u+0.5, 0.1.0*v, u w pm3d,\
D*u+1.5, 1.0*v, u w pm3d,\
D*u+2.5, 2.0*v, u w pm3d,\
D*u+3.5, 1.2*v, u w pm3d,\
D*u+4.5, 1.6*v, u w pm3d
unset ytics; unset ylabel
set palette model HSV function 0.75, 1.0-f(gray, A, B, C), (1.0+2.0*f(gray,A, B, C))/3.0
splot D*u+0.5,1.0+2.0*v,u w pm3d,\
D*u+1.5, 1.0+1.1*v, u w pm3d,\
D*u+2.5, 2.0+1.2*v, u w pm3d,\
D*u+3.5, 1.2+2.4*v, u w pm3d,\
D*u+4.5, 1.6+2.1*v, u w pm3d
set label 1 "First" at 0.5, -1.32 rotate by 45
set label 2 "Second" at 1.5, -1.32 rotate by 45
set label 3 "Third" at 2.5, -1.32 rotate by 45
set label 4 "Fourth" at 3.5, -1.32 rotate by 45
set label 5 "Fifth" at 4.5, -1.32 rotate by 45
set palette model HSV function 0.5, 1.0-f(gray, A, B, C), (1.0+2.0*f(gray,A, B, C))/3.0
splot D*u+0.5, 3.0+1.0*v, u w pm3d,\
D*u+1.5, 2.1+3.2*v, u w pm3d,\
D*u+2.5, 3.2+3.4*v, u w pm3d,\
D*u+3.5, 3.6+1.1*v, u w pm3d,\
D*u+4.5, 3.7+1.5*v, u w pm3d
unset multiplot

Most of the script has been analysed in the previous post, so I will point out the differences only. If you recall, when we plotted the bargraphs side by side, we had to shift the 'x' position, so that they didn't overlap. Now, since we want to stack them, we have to shift the 'y' position, and in particular, we have to shift it by a value that renders the bars just on top of each other. So, if you look at the three lines highlighted in blue, you will notice that the multiplier of the dummy variable 'v' is always added to the shift of 'v' in the consecutive plot, i.e., in the first plot, the multiplier of 'v' is 1 and the shift is 0, in the second plot, the multiplier of 'v' is 1.1 and the shift is 1, and in the third plot, the multiplier of 'v' is 3.2, and the shift is 2.1, which happens to be the sum of 1 and 1.1. These were the appropriate shift for the second column in the plot. In the rest of the script, we just have to follow the same recipe, and indeed, this is what happens. (I won't demonstrate this, though:-)

Obviously, one would like to do this automatically, and this can be done with the following gawk script:


#!/bin/bash

gawk  'BEGIN {i=0; max=0}
{
if($0!~/#/) {
label[i] = $1
sum=0
for(j=2;j<=NF;j++) {
      v[i,j-1] = $j
      sum+=$j
}
if(max<sum) max=sum
i++
}
}
END {
print "reset"
print "f(x, a, b, c) = a*exp(-(x-b)*(x-b)/c/c)"
print "A = 0.6; B = 0.5; C = 0.2; D=0.5"
print "unset key; unset colorbox; unset xtics; set pm3d map"
printf "set yrange [0:%f]\n", max
printf "set parametric; set urange [0:1]; set vrange [0:1]\n"
printf "set xrange [0:%f]\n", i+0.5
print "set multiplot"
print "set isosamples 100, 100"
print "set palette model HSV functions 2.0/3.0, 0.4-0.4*gray, (4.0+2.0*gray)/6.0"
printf "splot %f*u, %f*v, 2*u+v w pm3d\n", i+0.5, max
print "set isosamples 2, 2"
for(l=1;l<j-1; l++) {
 if(l==j-2) {
for(k=0;k<i;k++) {
      printf "set label %d \"%s\" at %f, -%f rotate by 45\n", k+1, label[k], k+0.5, max/5.0
}
 }
 if(l==2) print "unset ytics; unset ylabel"
 printf "set palette model HSV function %f, 1.0-f(gray, A, B, C), (1.0+2.0*f(gray, A, B, C))/3.0\n", 1.0-(l-1)/(j-1)
 printf "splot "
 for(k=0;k<i-1;k++) {
      printf "D*u+%f,%f+%f*v,u w pm3d,\\\n", k+0.5, base[k], v[k,l]
      base[k]+=v[k,l]
 }
 printf "D*u+%f,%f+%f*v,u w pm3d\n", i-0.5, base[i-1], v[i-1,l]
 base[i-1]+=v[i-1,l]
}
print "unset multiplot"

}' $1

It is very easy to turn this script into one that puts out normalised stacks. All we have to do is to divide the values by the appropriate sum, and re-set the 'y' range to [0:1.1] or something similar, shift the 'y' position of the labels (remember, it was defined as a function of 'max' above, but max is just 1.1 in the present case), and there is nothing else. This difference is highlighted in blue.


#!/bin/bash

gawk  'BEGIN {i=0; max=0}
{
if($0!~/#/) {
label[i] = $1
sum=0
for(j=2;j<=NF;j++) {
      v[i,j-1] = $j
      sum+=$j
}
if(sum!=0)
   for(j=2;j<=NF;j++) {
         v[i,j-1] /= sum
   }
i++
}
}
END {
print "reset"
print "f(x, a, b, c) = a*exp(-(x-b)*(x-b)/c/c)"
print "A = 0.6; B = 0.5; C = 0.2; D=0.5"
print "unset key; unset colorbox; unset xtics; set pm3d map"
print "set yrange [0:1.1]"
print "set parametric; set urange [0:1]; set vrange [0:1]"
printf "set xrange [0:%f]\n", i+0.5
print "set multiplot"
print "set isosamples 100, 100"
print "set palette model HSV functions 2.0/3.0, 0.4-0.4*gray, (4.0+2.0*gray)/6.0"
printf "splot %f*u, 1.1*v, 2*u+v w pm3d\n", i+0.5
print "set isosamples 2, 2"
for(l=1;l<j-1; l++) {
 if(l==j-2) {
for(k=0;k<i;k++) {
      printf "set label %d \"%s\" at %f, -0.2 rotate by 45\n", k+1, label[k], k+0.5
}
 }
 if(l==2) print "unset ytics; unset ylabel"
 printf "set palette model HSV function %f, 1.0-f(gray, A, B, C), (1.0+2.0*f(gray, A, B, C))/3.0\n", 1.0-(l-1)/(j-1)
 printf "splot "
 for(k=0;k<i-1;k++) {
      printf "D*u+%f,%f+%f*v,u w pm3d,\\\n", k+0.5, base[k], v[k,l]
      base[k]+=v[k,l]
 }
 printf "D*u+%f,%f+%f*v,u w pm3d\n", i-0.5, base[i-1], v[i-1,l]
 base[i-1]+=v[i-1,l]
}
print "unset multiplot"

}' $1

False 3D bargraphs in gnuplot

2009-06-14T05:06:00.000-07:00

Yesterday, we saw how cuboids representing bargraphs can easily be drawn in gnuplot. That implementation was completely 3D, so to speak. Today we will try something simpler. The idea is that in order to suggest that something is in 3D, it is not necessary to draw its 2D projection accurately, it is enough, if we can give the impression that the object is a three-dimensional one. Perhaps, the simplest way to do this is to colour the surface properly: by making some parts shiny, we can pretend that light is reflected from a curved surface. This is the trick that we are going to use today. First, we will take an easy example, and then later we can expand it. We will draw the following graph in a couple of lines:

As usual, we give the gnu script first, and then dissect it, step by step.

reset

f(x, a, b, c) = a*exp(-(x-b)*(x-b)/c/c)
A = 0.6
B = 0.5
C = 0.2

unset key
unset colorbox
unset xtics
set ytics 0,1,4
set parametric
set urange [0:0.5]
set vrange [0:1]

set pm3d map
set xrange [0:10]
set yrange [0:4]

set multiplot
set isosamples 100, 100
set palette model HSV function 2.0/3.0, 0.4-0.4*gray, (4.0+2.0*gray)/6.0
splot 20*u, 4*v, 2*u+v w pm3d

set isosamples 100, 2
set palette model HSV function 0.95, 1-f(gray, A, B, C), (1+2*f(gray, A, B, C))/3
splot u+1,v,u w pm3d,\
u+2,2*v,u w pm3d, \
u+3,3*v,u w pm3d, \
u+4,2*v,u w pm3d, \
u+5,3*v,u w pm3d, \
u+6,3.5*v,u w pm3d, \
u+7,1.4*v,u w pm3d, \
u+8,3*v,u w pm3d, \
u+9,1*v,u w pm3d

unset multiplot

First, we define the function according to which we are going to colour the bars. This is done by defining f(x,a,b,c) and the parameters A, B, and C. The the following couple of lines define the graph properties. The first relevant line is in blue, where we draw the background of our graph. We define the palette for pm3d using the HSV colour space. As in RGB, the colour is given as a triplet, the first member of which is the hue (in this case it will be something close to blue), the second is the saturation, which we defined to be whiter, if the value of gray is higher, and finally, the value. We then plot '20*u, 4*v, 2*u+v' on our map. By plotting the value 2*u+v, we will get a gradient that runs from bottom to top, left to right. You can modify this, if you want to change the direction in which the background becomes whiter.

Having drawn the background, we can plot the bars. This is done in the next step, where we always plot something like 'u+shift,value*v,u', where shift is the position of the bar, and value is its height. Note that before calling splot, we re-defined our palette, using the function f(x,a,b,c) from the beginning of the script. 'a' will set how white the colour becomes, 'b' is the position of the highest saturation value (or the centre of the pattern), and 'c' determines how tight the colour gradient is.

Now, let us see how this procedure can be automated. As previously, we will use a gawk script, which processes a data file with N+1 columns, the first setting the label, and the last N containing N values (i.e., instead of on set of bars, we will have N sets bars). Here is the script,


#!/bin/bash

gawk  'BEGIN {i=0; max=0}
 {
  if($0!~/#/) {
   label[i] = $1
   for(j=2;j<=NF;j++) {
         v[i,j-1] = $j
         if(max<v[i,j-1]) max=v[i,j-1]
   }
   i++
 }
 }
 END {
  print "reset"
  print "f(x, a, b, c) = a*exp(-(x-b)*(x-b)/c/c)"
  print "A = 0.6; B = 0.5; C = 0.2"
  print "unset key; unset colorbox; unset xtics; set pm3d map"
  printf "set yrange [0:%f]\n", max
  printf "set parametric; set urange [0:%f]; set vrange [0:1]\n", 1.0/(j-1)
  printf "set xrange [0:%d]\n", i+1
  print "set multiplot"
  print "set isosamples 100, 100"
  print "set palette model HSV functions 2.0/3.0, 0.4-0.4*gray, (4.0+2.0*gray)/6.0"
  printf "splot %d*u, %f*v, 2*u+v w pm3d\n", (i+1)*(j-1), max
  print "set isosamples 100, 2"
  for(l=1;l<j-1; l++) {
    if(l==j-2) {
   for(k=0;k<i;k++) {
         printf "set label %d \"%s\" at %f, -%f rotate by 45\n", k+1, label[k], k+0.5, max/5.0
  }
    }
    if(l>1) print "unset ytics; unset ylabel"
    printf "set palette model HSV function %f, 1.0-f(gray, A, B, C), (1.0+2.0*f(gray, A, B, C))/3.0\n", 1.0/l
    printf "splot "
    for(k=0;k<i;k++) {
         printf "u+%f,%f*v,u w pm3d,\\\n", k+l/(j-1), v[k,l]
    }
    printf "u+%f,%f*v,u w pm3d\n", i+l/j, v[i,l]
  }
  print "unset multiplot"
  
 }' $1

and here is the data file that I used to produce the figure below:

First 1.0 2.0 1.0
Second 1.0 1.1 3.2
Third 2.0 1.2 3.4
Fourth 1.2 2.4 1.1
Fifth 1.6 2.1 1.5

I believe, it should be fairly easy to hack the script to suit your needs. The points of interest are the colour of the background and the bars, and the position of the labels. Small tweaks in the script will give you just anything you would want. Next time I will come back to this script, and we will see how this can be modified to produce stacked bargraphs. So long!

3D bargraphs in gnuplot, one more time

2009-06-13T04:39:00.001-07:00

Some time ago, on 25th May, we discussed how we can draw bargraphs based on cylinders. What was good about them is that we gave the graph a 3D look by adding phong to the surface. The downside is that cylinders are not the conventional shapes to draw something like this. (As a matter of fact, there is one exception: when one wants to draw a bargraph in 2D, but give it a phong on the 2D surface, which gives the impression that what we are looking at is the projection of a cylinder. I will come to this point on a later date.) Instead, people plot a rectangular cuboid. Today we will try our hands at that. At the end of the day, we will have the following graph:

First, let us see how we create the figure above. The gnu script for drawing one cuboid is as follows:

set multiplot
set palette defined (0 1 0.5 0.5, 1 1 0.5 0.5)
splot u, v, B w pm3d

set palette defined (0 1 0.3 0.3, 1 1 0.3 0.3)
splot A+u, 0, v w pm3d

set palette defined (0 0.8 0 0, 1 0.8 0 0)
splot A+A, u, B*v w pm3d

unset multiplot

First, we define our palette. This should have only one colour per face, because we do not want to add any gradient, phong etc. (If you insist, you can add that very easily.) Then we plot the top of the cuboid, at a z-level of B, and shifted horizontally by A. In the next step, we re-define our palette, which will now have a slightly darker shade of red, and plot the front surface. In the final step, we move to the third visible face, make it darker, and plot it. In this way, we have drawn one cuboid. The rest is nothing but the repetition of these three steps, and the general set-up of the graph. Obviously, by changing the order of the palette definitions, one can change the direction of the lighting.

Instead of copying the complete gnu script here, I will just give the gawk script, which processes a file of 2 columns, the first being the label of the bars, and the second containing the actual values. As usual, you can call the script from gnuplot without writing the file to disc as

load "< gawk-script.sh my_data_file.dat"

So, this is what we have:

#!/bin/bash

gawk  'BEGIN {i=0; max=0}
 {
  if($0!~/#/) {
   label[i] = $1
   v[i] = $2
   if(max<v[i]) max=v[i];
   i++
 }
 }
 END {
  print "reset"
  print "set view 60, 20; set parametric; set isosample 2, 2"
  print "unset key; unset colorbox"
  print "set ticslevel 0"
  print "set urange [0:0.5]; set vrange [0:1]"
  printf "set xrange [0:%d]; set yrange [-%f:%f]; set zrange [0:%f]\n", i, i/2, i/2, max
  print "set multiplot"
  print "set border 1+2+4+8+16+32+64+256+512; unset ytics; unset xtics"
  print "set palette model RGB functions 0.9, 0.9,0.95"
  printf "splot %d*u, -%f+%d*v, 0 w pm3d\n", 2*i, i/2, i
  print "unset border; unset xtics; unset ytics; unset ztics"
  print "set palette model RGB functions 1, 254.0/255.0, 189.0/255.0"
  printf "splot 0, -%f+%d*v, %f*u w pm3d, %d*u, %f, %f*v w pm3d\n", i/2, i, 2*max, 2*i, i/2, max

  print "set palette defined (0 1 0.5 0.5, 1 1 0.5 0.5)"
  printf "splot "
  for(j=0;j<i-1;j++) {
   printf "%d+u, v, %f w pm3d,\\\n", j, v[j]
  }
  if(i>1) printf "%d+u, v, %f w pm3d\n", i-1, v[i-1]
  
  print "set palette defined (0 1 0.3 0.3, 1 1 0.3 0.3)"
  printf "splot "
  for(j=0;j<i-1;j++) {
   printf "%d+0.5, 2*u, v*%f w pm3d,\\\n", j, v[j]
  }
  if(i>1) printf "%d+0.5, 2*u, v*%f w pm3d\n", i-1, v[i-1]
  
  print "set palette defined (0 0.8 0 0, 1 0.8 0 0)"
  for(j=0;j<i;j++) {
   printf "set label %d \"%s\" at %f, %f rotate by 40\n", j+1, label[j], j+0.3, -i/2+1
  }
  printf "splot "
  for(j=0;j<i-1;j++) {
   printf "%d+u, 0, v*%f w pm3d,\\\n", j, v[j]
  }
  
  if(i>1) printf "%d+u, 0, v*%f w pm3d\n", i-1, v[i-1]
 
  print "unset multiplot"
  
 }' $1

The generalisation to more than 1 data column should be trivial: you have only got to choose your colours, and re-size the bars in the x-y directions accordingly.

Wall charts with gnuplot

2009-06-07T04:42:00.000-07:00

Today we will discuss how to make a wall chart with gnuplot. There is a simpler version of this, the so-called fence plot, that you can find on gnuplot's surface demo page, and we will build on this. By the end of our adventure, we will have produced this graph:

So, let us see what is to be done. First, I will just give the gnu script, and then discuss segments of it.

reset

A = 4
B = 3
C = 1
unset key
unset colorbox
f(x,y) = abs(sin(x*x+y*y)/(x*x+y*y)) # We plot this function
set isosample 2, 30
set parametric
set urange [-B:B]
set vrange [0:1]
set xrange [-A:A]
set yrange [-B:B]
set zrange [0:C]
set ticslevel 0 # We need this, in order to shift the 0 of the z axis

dx = 0.3 # Thickness of the walls

set multiplot
set border 1+2+4+8+16+32+64+256+512
# Plot the "background" first
set palette model RGB functions 0.9, 0.9,0.95
splot -A+2*A*v, B*u,0 w pm3d

unset border; unset xtics; unset ytics; unset ztics
set palette model RGB functions 1, 254.0/255.0, 189.0/255.0
splot -A, B*u, C*v w pm3d, A*u, B, C*v w pm3d

x0 = -A+1
y0 = 0.0
set palette model RGB functions 0, gray/2+0.5, 0.0
splot x0-dx, u, f(u,y0) w l lt -1 lw 0.2, \
x0-dx*v, u, f(u,y0) w pm3d,\
x0, u, f(u,y0)*v with pm3d, \
x0, u, f(u,y0) w l lt -1 lw 0.2

x0 = x0+1
y0 = y0+0.5
set palette model RGB functions gray/2+0.5, 0.0, 0.0
splot x0-dx, u, f(u,y0) w l lt -1 lw 0.2, \
x0-dx*v, u, f(u,y0) w pm3d,\
x0, u, f(u,y0)*v with pm3d, \
x0, u, f(u,y0) w l lt -1 lw 0.2

x0 = x0+1
y0 = y0+0.5
set palette model RGB functions 0, 0, gray/2+0.5
splot x0-dx, u, f(u,y0) w l lt -1 lw 0.2, \
x0-dx*v, u, f(u,y0) w pm3d,\
x0, u, f(u,y0)*v with pm3d, \
x0, u, f(u,y0) w l lt -1 lw 0.2

x0 = x0+1
y0 = y0+0.5
set palette model RGB functions gray/2+0.5, gray/2+0.5, 0
splot x0-dx, u, f(u,y0) w l lt -1 lw 0.2, \
x0-dx*v, u, f(u,y0) w pm3d,\
x0, u, f(u,y0)*v with pm3d, \
x0, u, f(u,y0) w l lt -1 lw 0.2

x0 = x0+1
y0 = y0+0.5
set palette model RGB functions gray/2+0.5, 0, gray/2+0.5
splot x0-dx, u, f(u,y0) w l lt -1 lw 0.2, \
x0-dx*v, u, f(u,y0) w pm3d,\
x0, u, f(u,y0)*v with pm3d, \
x0, u, f(u,y0) w l lt -1 lw 0.2

x0 = x0+1
y0 = y0+0.5
set palette model RGB functions 0, gray/2+0.5, gray/2+0.5
splot x0-dx, u, f(u,y0) w l lt -1 lw 0.2, \
x0-dx*v, u, f(u,y0) w pm3d,\
x0, u, f(u,y0)*v with pm3d, \
x0, u, f(u,y0) w l lt -1 lw 0.2

x0 = x0+1
y0 = y0+0.5
set palette model RGB functions gray/3+0.3, 0, 0
splot x0-dx, u, f(u,y0) w l lt -1 lw 0.2, \
x0-dx*v, u, f(u,y0) w pm3d,\
x0, u, f(u,y0)*v with pm3d, \
x0, u, f(u,y0) w l lt -1 lw 0.2
unset multiplot

In the first part, we set up the plotting ranges and the function that we want to plot. Note that in order to shift the 0 of the z axis, we have explicitly got to call 'set ticslevel 0'. We can also notice that in the graph above, the top surface is invariant under translations along the x axis, therefore, we can reduce the sampling frequency in that direction. This shouldn't matter, if you make a bitmap file, but it is going to save you a lot of space, if you produce a vector output, e.g., eps or svg.

Having finished the set-up, we plot the background, i.e., we will just plot three planes, the x-y (this one in gray), the y-z, and the x-z. These latter ones are in yellow. We set the colours in the colour scheme of the appropriate plots, by calling 'set palette model RGB functions 0.9, 0.9,0.95', which produces a shade of gray with a bluish tinge. Also note that we unset the border, and the tics immediately after the first plot. The reason for this is that in this way, we can avoid having gnuplot plot the axes and the tics more than once. By the way, before the first plot, we set the border as 'set border 1+2+4+8+16+32+64+256+512', which draws the frame along all, but the three front edges. You can learn how to set the border by reading the output of the command

?border

After having plotted the 'background' of the graph, we plot the first wall. In order to give it some 3D lookout, we have got to plot actually 4 different things: the vertical wall facing us, the top of the body, and the two relevant edges. The first two plots are surfaces, while the latter two are lines, so we call appropriate modifiers to splot, namely, 'with pm3d' and 'with lines'. We have also got to specify the colour of the walls, which we did by issuing 'set palette model RGB functions 0, gray/2+0.5, 0.0' before the first plot. The rest of the script is nothing but the repetition of these seven lines. Letting x0 = x0 + 1, we move the walls along the x axis, while letting y0 = y0 + 0.5, we plot a new function each time. Since we want each wall to have a different colour, we also need to change the palette in each step.

There are two remarks that we can make. The first is that the order of plotting the wall, the top and the two edges depends on the view, so don't be surprised, if you rotate the figure and all of a sudden, it looks messy! The second comment is that you can easily 'close' the graphs, by adding one parametric plot to each of the walls, which would be something like this:


x0-dx*v, -B, f(x0,-B)*(B+u)/(2*B) w pm3d

This is it for today!

A comment on phonged surfaces in Gnuplot - how can the colour space be expanded?

2009-06-04T14:41:00.000-07:00

A couple of days ago (on the 24th of May, to be more specific), we discussed a way to add phong to a plot. However, that method worked only for a single colour (red in that particular case), but it would not do any good for a palette. Palette is the name of the colour space that is used in pm3d to colour surfaces. At the end of the post, I also showed a phonged surface with a general surface, and I must say that that wasn't very appealing. The problem is, if you take a look at this graph

you will immediately notice that apart from the phong, not too much can be seen, for the simple reason that the surface is not coloured, so it is very hard to get a feeling of the heights and depths. Using pm3d, we can colour the surface, but we can't add phong. So, the question is, how do we combine the figure above with this one:

While it appears to be trivial, it turns out to be rather hard. But will manage, don't worry, and at the end of the day, we will have figures that look like these:

Understanding the colouring scheme

To begin, we have got to understand how the colouring in pm3d is done. When you plot without specifying the colour scheme, you are given the default, which looks like this (look at the colour box at the bottom, don't worry about the curves)

If you dig a bit deeper, e.g., by issuing the help command


?rgbformulae

you will be presented with cryptic messages, like


Some nice schemes in RGB color space
7,5,15   ... traditional pm3d (black-blue-red-yellow)
3,11,6   ... green-red-violet
23,28,3  ... ocean (green-blue-white); try also all other permutations
21,22,23 ... hot (black-red-yellow-white)
30,31,32 ... color printable on gray (black-blue-violet-yellow-white)
33,13,10 ... rainbow (blue-green-yellow-red)
34,35,36 ... AFM hot (black-red-yellow-white)

So, what does this all mean?

When invoking pm3d, first, the colour range (cbrange) is calculated (or specified), and then, using the colour range, all values of the function z(x,y) are mapped to the [0:1] interval. Using this value, which is called gray, the three components of the colour palette are calculated through three functions, red(gray), green(gray) and blue(gray). When you invoke the command


set palette rgbformulae 7, 5, 15

(This is the default, so you actually haven't got to invoke this command...), you assign function 7 to red(gray), function 5 to green(gray) and function 15 to blue(gray). But then, how on Earth does one figure out which function to use? Well, the easy solution is to issue


show palette rgbformulae

upon which, this list is shown:


* there are 37 available rgb color mapping formulae:
0: 0               1: 0.5             2: 1
3: x               4: x^2             5: x^3
6: x^4             7: sqrt(x)         8: sqrt(sqrt(x))
9: sin(90x)       10: cos(90x)       11: |x-0.5|
12: (2x-1)^2       13: sin(180x)      14: |cos(180x)|
15: sin(360x)      16: cos(360x)      17: |sin(360x)|
18: |cos(360x)|    19: |sin(720x)|    20: |cos(720x)|
21: 3x             22: 3x-1           23: 3x-2
24: |3x-1|         25: |3x-2|         26: (3x-1)/2
27: (3x-2)/2       28: |(3x-1)/2|     29: |(3x-2)/2|
30: x/0.32-0.78125 31: 2*x-0.84       32: 4x;1;-2x+1.84;x/0.08-11.5
33: |2*x - 0.5|    34: 2*x            35: 2*x - 0.5
36: 2*x - 1
* negative numbers mean inverted=negative colour component
* thus the ranges in `set pm3d rgbformulae' are -36..36

So, with these in mind, the default pm3d colouring scheme is sqrt(gray) for red, x^3 for green and sin(360x) for blue. These three functions are shown in the previous figure.

So far, so good, but how do we turn all this into a phong on a general surface? If you recall from the post on phonged surfaces, what we did was to make the colour whiter when the value of the Gaussian function was high, and redder, when it was small. This worked very well for a single colour, because in that case, we simply fixed the red function at 1, and changed the value of the green and blue components according to the value of the Gaussian. But if we want to retain the pm3d colouring and apply phong at the same time, we have a bit of a problem here, because all of a sudden, the colour of a point on the surface will depend on two quantities: the value if z(x,y), and the value of f(x,y,z), which is our Gaussian function. However, the colour functions take only one argument, and nothing more. How will we parametrize a 2D space with only one variable? Well, we will discretise it in a smart way. Let us suppose, that we want to have only 3 white levels, and let us suppose, for a moment, that we are not restricted to the [0:1] interval for gray, but it can run between 0 and 4. Then, using only the red function for now, we will take the following scheme:


red(x) = sqrt(x)*3/3 + 0, if x is in [0:1]
red(x) = sqrt(x-1)*2/3 + 1/3, if x is in [1:2]
red(x) = sqrt(x-2)/3 + 2/3, if x is in [2:3]
red(x) = sqrt(x-3)*0 + 3/3, if x is in [3:4]

So, what happens here is that the red component has the same functional behaviour (namely, sqrt(x)) in all four intervals, but the base value is shifted upwards (to the white), while the amplitude is decreased in such a way that red(x) is always in the [0:1] interval. Obviously, we can play the same trick with the green and blue components. What we should notice is that in each interval, the argument of the original function is the fractional part of x, and that the coefficient of sqrt(x) plus the base line adds up to one. Therefore, we can conveniently write


red(x) = RED(frac(x))*(1-floor(x)/N) + floor(x)/N

where frac(x) is just the fractional part of x,


frac(x) = x - floor(x)

and RED(x) is the original colouring function (sqrt(x) in the example above). N was the number of white levels (3 above) that we wanted to have. Now we have almost everything at our disposal in order to add phong to an arbitrary surface. The last thing to do is to make sure that the argument of red(x) etc. falls into the [0:1] interval. We achieve this buy expanding the colour range with respect to the zrange (which is the interval of possible z values). With these in mind, the gnu file would read as follows:


reset

unset colorbox
unset key
unset border
unset tics
f(x,y) = sin(x*x+y*y)/(x*x+y*y)

A = 100.0
zi = -0.4
za = 1.0
ci = zi
ca = zi + (za-zi)*A

set isosample 150, 150
set cbrange [ci:ca]
set xrange [-3:3]
set yrange [-3:3]
set zrange [zi:za]

set table 'phong2.dat'
splot f(x,y)
unset table

# These are the functions of the original colour palette
r(x) = sqrt(x)
g(x) = x**3.0
b(x) = q(x<0.5?sin(x*2.0*pi):0.0)

# The Gaussian phong
h(x,a,b) = exp(-(x-a)*(x-a)/b/b)

# These are helper functions
R(x) = A*x - A*floor(x) # Fractional part of x
L(x) = floor(x)         # This is just a short-hand for the floor function
F(x) = L(A*x)/(ca-ci)   # The floor function divided by the colour range

# This is where we want to place the white spot
x0 = -0.3
y0 = -1.0
z0 = f(x0, y0)

# And these are the new colour scheme
red(gray) = r(R(gray))*(1.0-F(gray))+F(gray)
green(gray) = green(R(gray))*(1.0-F(gray))+F(gray)
blue(gray) = blue(R(gray))*(1.0-F(gray))+F(gray)

set palette model RGB functions red(gray), green(grey), blue(gray)
splot 'phong2.dat' u 1:2:3:((A-1)*(za-zi)*h($3,z0,0.08)*h($1,x0,0.3)*h($2,y0,0.3)+$3)) w pm3d

That is, first we decide that we want to have 100 white levels, this is A = 100. Then we fix the zrange according to the function that we want to plot, and after that, we specify the colour range. Here we shouldn't forget that have got to stretch the colour range to 100 times the zrange. Next, we identify the functions that we will need. For a start, we begin with the default pm3d scheme, whose functions we have already discussed above. Then a couple of helper functions follows, with the definition of the new colour scheme. Immediately before plotting the surface, we set the palette, and then we call the file that we produced a couple of lines earlier. We plot the surface using the x, y, z values, and for the colouring, we stretch the range of gray as we discussed above: the value of z(x,y) (this is the 3rd column, $3) will determine the base colour, and the value of the Gaussian will determine how much white tinge this base colour will have.

Well, there is not too much more to it. Once you decided your colour scheme and the number of white levels, the rest should be plain sailing.