Walking Randomly

MATLAB contains a function called pdist that calculates the ‘Pairwise distance between pairs of objects’. Typical usage is

X=rand(10,2);
dists=pdist(X,'euclidean');

It’s a nice function but the problem with it is that it is part of the Statistics Toolbox and that costs extra. I was recently approached by a user who needed access to the pdist function but all of the statistics toolbox license tokens on our network were in use and this led to the error message

??? License checkout failed.
License Manager Error -4
Maximum number of users for Statistics_Toolbox reached.
Try again later.
To see a list of current users use the lmstat utility or contact your License Administrator

One option, of course, is to buy more licenses for the statistics toolbox but there is another way. You may have heard of GNU Octave, a free,open-source MATLAB-like program that has been in development for many years.  Well, Octave has a sister project called Octave-Forge which aims to provide a set of free toolboxes for Octave.  It turns out that not only does Octave-forge contain a statistics toolbox but that toolbox contains an pdist function.  I wondered how hard it would be to take Octave-forge’s pdist function and modify it so that it ran on MATLAB.

Not very!  There is a script called oct2mat that is designed to automate some of the translation but I chose not to use it – I prefer to get my hands dirty you see.  I named the resulting function octave_pdist to help clearly identify the fact that you are using an Octave function rather than a  MATLAB function.  This may matter if one or the other turns out to have bugs.  It appears to work rather nicely:

dists_oct = octave_pdist(X,'euclidean');
% let's see if it agrees with the stats toolbox version
all( abs(dists_oct-dists)<1e-10)

ans =
     1

Let’s look at timings on a slightly bigger problem.

>> X=rand(1000,2);
>> tic;matdists=pdist(X,'euclidean');toc
Elapsed time is 0.018972 seconds.
>> tic;octdists=octave_pdist(X,'euclidean');toc
Elapsed time is 6.644317 seconds.

Uh-oh! The Octave version is 350 times slower (for this problem) than the MATLAB version. Ouch! As far as I can tell, this isn’t down to my dodgy porting efforts, the original Octave pdist really does take that long on my machine (Ubuntu 9.10, Octave 3.0.5).

This was far too slow to be of practical use and we didn’t want to be modifying algorithms so we ditched this function and went with the NAG Toolbox for MATLAB instead (routine g03ea in case you are interested) since Manchester effectively has an infinite number of licenses for that product.

If,however, you’d like to play with my MATLAB port of Octave’s pdist then download it below.

• octave_pdist.m makes use of some functions in the excellent NaN Toolbox so you will need to download and install that package first.

Did you make any New Year’s resolutions this year?  If you did then who will they help?  Just you?  Your family?  Your students?  The whole world?  If I am being honest then I have to say that most of the new-year’s resolutions I have made over the years tend to focus on myself because at my very core I am a bit selfish.  So, my resolutions tend to be things like “I want to get fitter”, “I’ll not stay late at work so much” or “I want to learn more Python programming.”

If I keep these resolutions then I’m going to be fitter, more knowledgeable and have a better work-life balance.  So far so selfish!

Over the last few days though I have made a rather different sort of new-years resolution.  Yes, I admit that it’s a bit late but why limit change to an essentially arbitrary date?  My new new year’s resolution is to give a little back to the community that has given me so much – the community of organisations and individuals who provide me with software – either for free or for such a trivial amount of money that it may as well be free.

Now I am not a rich man so I can’t give away great wads of cash and although I am a programmer I have neither the time nor the knowledge to provide significant amounts of code to any particualr project.  So what can I do?

Donations

Well, although I am not minted, I can easily afford the occasional small donation or two so I will start making them.  It’s already begun with my Sage bounty hunt and will continue with other projects over the year.  Another small donation I made recently was to buy the ‘Premium’ version of Aldiko – a great ebook reader for Android smartphones.  Aldiko is a free piece of software, has been downloaded by tens of thousands of users and is used by me on an almost daily basis.  I noticed that they had a ‘Premium’ version available for $1.99 but it turns out that it is identical to the free version.  The $1.99 simply represents a donation to the developers and it’s a donation I made without hesitation.  Doing the right thing for less than two dollars – getting the warm fuzzies has never been easier!

Bug reports and feedback

Thanks to my job and to running Walking Randomly I get to see how researchers, teachers, students and average joes use mathematical software quite a lot.  I get told about bugs, about feature wish lists, about gripes with licensing, performance issues…the list goes on and on.  The best way to get bugs fixed is to report them – first tell the developers and, if the bug is interesting/severe enough, tell the world.  I do this plenty with commercial software but I am going to make the effort more with free software from now on.  Feedback is part of the lifeblood of free software and developers need both the good and the bad.

Did you try out Octave, Maxima, Sage or Smath Studio and it didn’t work out?  Why didn’t it?  What did these packages have missing that forced you to turn to alternatives?  Try to be specific; saying ‘I tried it and it sucks’ is a rubbish piece of feedback because it’s just an opinion and gives the developers nothing to do.  Saying something like ‘I tried to calculate XYZ and it gave the following incorrect result whereas it should have given ABC’ is MUCH more productive.

Tutorials, examples and documentation

One comment I have heard over and over again from people who have tried free mathematical software and then turned their back on it is that the ‘documentation isn’t good enough’.  These people want more tutorials, more examples, more explanations and just more and better documentation.  Do you like writing?  Do you like fiddling with math software?  I do and so I intend on giving as many examples and tutorials as I can.  I also have a (moderately) successful blog so I can provide a publishing outlet for people who want to write such things but don’t want to start a blog of their own.  This has already begun too with Greg Astley’s tutorial on how to plot direction fields for ODE’s in maxima.  Contact me if you are interested in doing this yourself and we’ll discuss it.

Talk

I like to talk.  Many people who know me personally would probably go so far as to say I talk too much but I can use this to help towards my new resolution too.  I’m going to give short talks, demonstrations and seminars on free mathematical software to interested people over the coming year via various fora.  Maybe you could too?

So, in summary, I plan to do the following to give back to free software over the next decade and I invite you to do the same.

• Give small, direct donations to some of my favourite open source and free software projects
• Set up bounty hunts for particular features I want in various packages
• Buy donation versions of Android software whenever possible
• Publish as many examples of using software such as Sage, octave and maxima as I can
• Help write tutorials and documentation
• Give talks to help spread the word

I’ll admit that none of this will change the world but it will possibly help a few more people than “I resolve to get myself fitter.”

This is the first post on Walking Randomly that isn’t written by me!  I’ve been thinking about having guest writers on here for quite some time and when I first saw the tutorial below (written by Manchester University undergraduate, Gregory Astley) I knew that the time had finally come.  Greg is a student in Professor Bill Lionheart’s course entitled Asymptotic Expansions and Perturbation Methods where Mathematica is the software of choice.

Now, students at Manchester can use Mathematica on thousands of machines all over campus but we do not offer it for use on their personal machines.  So, when Greg decided to write up his lecture notes in Latex he needed to come up with an alternative way of producing all of the plots and he used the free, open source computer algebra system – Maxima.  I was very impressed with the quality of the plots that he produced and so I asked him if he would mind writing up a tutorial and he did so in fine style.  So, over to Greg….

This is a short tutorial on how to get up and running with the “plotdf” function for plotting direction fields/trajectories for 1st order autonomous ODEs in Maxima. My immediate disclaimer is that I am by no means an expert user!!! Furthermore, I apologise to those who have some experience with this program but I think the best way to proceed is to assume no prior knowledge of using this program or computer algebra systems in general.  Experienced users (surely more so than me) may feel free to skip the ‘boring’ parts.

Firstly, to those who may not know, this is a *free* (in both the “costs no money”, and “it’s open source” senses of the word) computer algebra system that is available to download on Windows, Mac OS, and Linux platforms from  http://maxima.sourceforge.net and it is well documented.

There are a number of different themes or GUIs that you can use with the program but I’ll assume we’re working with the somewhat basic “Xmaxima” shell.

Install and open up the program as you would do with any other and you will be greeted by the following screen.

You are meant to type in the top most window next to (%i1) (input 1)
We first have to load the plotdf package (it isn’t loaded by default) so type:

load("plotdf");

and then press return (don’t forget the semi-colon or else nothing will happen (apart from a new line!)). it should respond with something similar to:

(%i1)  load("plotdf");
(%o1)  /usr/share/maxima/5.17.1/share/dynamics/plotdf.lisp
(%i2)

We will now race ahead and do our first plot. Keeping things simple for now we’ll do a phase plane plot for dx/dt = y, dy/dt = x, type:

plotdf([y,x]);

you should see something like this:

This is the Openmath plot window, (there are other plotting environments like Gnuplot but this function works only with Openmath) Notice that my pointer is directly below a red trajectory. These plots are interactive, you can draw other trajectories by clicking elsewhere. Try this.  Hit the “replot” button and it will redraw the direction field with just your last trajectory.

Before exploring any other options I want to purposely type some bad input and show how to fix things when it throws an error or gets stuck. Type

plotdf(y,x);

it should return

(%i3) plotdf(x,y);
Maxima encountered a Lisp error:

Error in SECOND [or a callee]: $Y is not of type LIST.

Automatically continuing.
To enable the Lisp debugger set *debugger-hook* to nil.
(%i4)

We forgot to put our functions for dx/dt,dy/dt in a list (square brackets). This is a reasonably safe error in that it tells us it isn’t happy and lets us continue.
Now type

plotdf([x.y]);

you should see something similar to

The problem this time was that we typed a dot instead of a comma (easily done), but worryingly when we close this message box and the blank plot the program will not process any commands. This can be fixed by clicking on the following in the xmaxima menu

file >> interrupt

where after telling you it encountered an error it should allow you to continue. One more; type

plotdf([2y,x]);

It should return with

  (%i5) plotdf([2y,x]);
  Incorrect syntax: Y is not an infix operator
  plotdf([2y,
           ^
  (%i5)

This time we forgot to put a binary operation such as * or + between 2 and y. If you come up with any other errors and the interrupt command is of no use you can still partially salvage things via

  file >> restart

but you will, in this case, have to load plotdf again. (mercifully you can go to the line where you first typed it and press return (as with other commands you might have done))
I will now demonstrate some more "contrived" plots (for absolutely no purpose other than to shamelessly give a (very) small gallery of standard functions/operations etc... for the novice user) there is no need to type the last four unless you want to see what happens by changing constants/parameters, they're the same plot

  plotdf([2*y-%e^(3/2)+cos((%pi/2)*x),log(abs(x))-%i^2*y]);
  plotdf([integrate(2*y,y)/y,diff((1/2)*x^2,x)]);
  plotdf([(3/%pi)*acos(1/2)*y,(2/sqrt(%pi))*x*integrate(exp(-x^2),x,0,inf)]);
  plotdf([floor(1.43)*y,ceiling(.35)*x]);
  plotdf([imagpart(x+%i*y),(sum(x^n,n,0,2)-sum(x^j,j,0,1))/x]);

I could go on…notice that the constants pi, e, and i are preceded by “%”. This tells maxima that they are known constants as opposed to symbols you happened to call “pi”, “e”, and “i”. Also, noting that the default range for x and y is (-10,10) in both cases; feel free to replot the first of those five without wrapping x inside “abs()” (inside the logarithm that is). remember file >> interrupt afterwards!

Now I will introduce you to some more of the parameters you can plug into “plotdf(…)”. close any plot windows and type

plotdf([x,y],[x,1,5],[y,-5,5]);

You should notice that x now ranges from 1 to 5, whilst y ranges from -5 to 5. There is nothing special about these numbers, we could have chosen any *real* numbers we liked. You can also use different symbols for your variables instead of x or y. Try

plotdf([v,u],[u,v]);

Note that I’ve declared u and v as variables in the second list. I will now purposely do something wrong again. Assign the value 5 to x by typing

x:5;

then type

plotdf([y,x]);

This time maxima won’t throw an error because syntactically you haven’t done anything wrong, you merely told it to do

plotdf([y,5]);

as opposed to what you really wanted which is

plotdf([y,x]);

Surprisingly to me (discovered as I’m writing this), changing the names of your variables like we did above won’t save you since it seems to treat new symbols as merely placeholders for it’s favourite symbols x and y. To get round this type

kill(x);

and this will put x back to what it originally was (the symbol x as opposed to 5).
You don’t actually have to provide expressions for dx/dt and dy/dt, you might instead know dy/dx and you can generate phaseplots by typing say

plotdf([x/y]);

In this case we didn’t actually need the square brackets because we are providing only one parameter: dy/dx (x will be set to t by maxima giving dx/dt = 1, and dy/dt = dy/dx = x/y)

A number of parameters can be changed from within the openmath window. Type

plotdf([y,x],[trajectory_at,-1,-2],[direction,forward]);

and then go into config. The screen you get should look something like this:

from here you can change the expressions for dx/dt, dy/dt, dy/dx, you can change colours and directions of trajectories (choices of forward, backward, both), change colours for direction arrows and orthogonal curves, starting points for trajectories (the two numbers separated by a space here, not a comma), time increments for t, number of steps used to create an integral curve. You can also look at an integral plots for x(t) and y(t) corresponding to the starting point given (or clicked) by hitting the “plot vs t” button. You can also zoom in or out by hitting the “zoom” button and clicking (or shift+clicking to unzoom), doing something else like hitting the config button will allow you to quit being in zoom mode click for trajectories again. (there might be a better way of doing this btw) You can also save these plots as eps files (you can then tweak these in other vector graphics based programs like Inkscape (free) or Adobe Illustrator etc..)

Interactive sliders

There are many permutations of things you can do (and you will surely find some of these by having a play) but my particular favourite is the option to use sliders allowing you to vary a parameter interactively and seeing what happens to the trajectories without constant replotting. ie:

plotdf([a*y,b*x],[sliders,"a=-1:3,b=-1:3"]);

Hopefully, this has been enough to get people started, and for what it’s worth, the help file (though using xmaxima, you’ll find this in the web-browser version) for this topic has a pretty good gallery of different plots and other parameters I didn’t mention.

just to throw in one last thing in the spirit of experimentation, is the following set of commands:

A:matrix([0,1],[1,0]);
B:matrix([x],[y]);
plotdf([A[1].B,A[2].B);

which is another way of doing the same old

plotdf([y,x]);

where here I’ve made a 2×2 matrix A, a 2×1 matrix B, with A[1], A[2] denoting rows 1 and 2 of A respectively and matrix multiplied the rows of A by B (using A[1].B, A[2].B) to serve as inputs for dx/dt and dy/dt

Tutorial written by Gregory Astley, Mathematics Undergraduate at The University of Manchester.

It’s been a great week for Open Source maths software this week with new releases of both Maxima and Scilab.

Maxima is a computer algebra system written in lisp that runs on most operating systems including Linux, Mac OS X and Windows.  The latest version, Maxima 5.20.1,  was released on Sunday 13th December and the full list of changes can be found here.  Highlights include improvements to the calculation of special functions, faster fourier transform routines and a general mechanism for functions to distribute over operators.

This last item is pretty cool since it allows you to map functions over lists in a similar manner to Mathematica.  For example, if you apply the sin function to a list in Maxima then you’ll now get a list containing the sines of each element of that list:

(%i1) sin([1,2,3]);
(%o1) [sin(1),sin(2),sin(3)]

Here are some examples of this new functionality for two-argument functions such as expintegral_e (taken from a post of Dieter Kaiser’s) include :

(%i2) expintegral_e([1,2,3],x);
(%o2) [expintegral_e(1,x),expintegral_e(2,x),expintegral_e(3,x)]

(%i3) expintegral_e(1,[x,y,z]);
(%o3) [expintegral_e(1,x),expintegral_e(1,y),expintegral_e(1,z)]

(%i4) expintegral_e([1,2,3],[x,y,z]);
(%o4) [[expintegral_e(1,x),expintegral_e(1,y),expintegral_e(1,z)],
       [expintegral_e(2,x),expintegral_e(2,y),expintegral_e(2,z)],
       [expintegral_e(3,x),expintegral_e(3,y),expintegral_e(3,z)]]

Head over to sourceforge to download Maxima for your operating system of choice. Now, onto Scilab….

Scilab is a fantastic open source numerical mathematics environment that always seems to come top of the list when people start discussing free MATLAB alternatives.  December 17th saw the release of Scilab 5.20 and the list of changes is epic (warning:pdf file).

The first new Scilab feature that caught my eye is the addition of a module called Xcos.  Xcos 1.0 is based on Scicos 4.3 which is a “graphical dynamical system modeler and simulator.”  Now I’ve not used any of this stuff before but as far as I can tell it could be considered as a free version of Simulink.  A bit of googling turned up a paper by M.G.J.M. Maassen, a student at the Eindoven University of Technology, who looked at the issue of migrating from Simulink to Scicos with respect to real time programs.  Written in 2006, this paper is slightly out of date but is a great start for anyone who is considering moving away from Simulink.

Along with other free mathematical software applications such as Octave, Sage, Python/Numpy and freemat, these new releases demonstrate that the world of free, open source mathematical applications has never looked better.

One of my plans for next year is to start giving short talks and tutorials about scientific Python to various groups around the University of Manchester.  Ever since I attended EuroSciPy in Leipzig earlier this year I have been excited about the possibilities offered by Python for students, researchers and teachers.  I genuinely believe that Python and SAGE will do for MATLAB and Mathematica users what R has done for users of SPSS and Stata.  It’s the future…I’ve tasted it!  Here are a few links to others who agree with me.

• Scientific Python tools use rising in education – a blog post from Fernando Perez
• A hands-on 2-day workshop at UC Berkeley – includes videotaped lectures
• Python for High Performance Computing (HPC)
• Should I switch to Python – Are you a MATLAB user considering the switch?
• Python in the Scientific World – From Guido van Rossum – Pythons benevolent dictator

Earlier this week I mentioned my plans to a colleague of mine over coffee and he said “Cool, I didn’t realise that you were a Python expert!”

Of course I’m NOT an expert on Python but it turns out that you just don’t have to be an expert in Python in order to get useful stuff done and THAT is my point!

There have been a couple of blog posts recently that have focused on creating interactive demonstrations for the discrete logistic equation – a highly simplified model of population growth that is often used as an example of how chaotic solutions can arise in simple systems.  The first blog post over at Division by Zero used the free GeoGebra package to create the demonstration and a follow up post over at MathRecreation used a proprietary package called Fathom (which I have to confess I have never heard of – drop me a line if you have used it and like it).  Finally, there is also a Mathematica demonstration of the discrete logistic equation over at the Wolfram Demonstrations project.

I figured that one more demonstration wouldn’t hurt so I coded it up in SAGE – a free open source mathematical package that has a level of power on par with Mathematica or MATLAB.  Here’s the code (click here if you’d prefer to download it as a file).

def newpop(m,prevpop):
    return m*prevpop*(1-prevpop)

def populationhistory(startpop,m,length):
    history = [startpop]
    for i in range(length):
        history.append( newpop(m,history[i]) )
    return history

@interact
def _( m=slider(0.05,5,0.05,default=1.75,label='Malthus Factor') ):
    myplot=list_plot( populationhistory(0.1,m,20) ,plotjoined=True,marker='o',ymin=0,ymax=1)
    myplot.show()

Here’s a screenshot for a Malthus Factor of 1.75

and here’s one for a Malthus Factor of 3.1

There are now so many different ways to easily make interactive mathematical demonstrations that there really is no excuse not to use them.

Update (11th December 2009)

As Harald Schilly points out, you can make a nice fractal out of this by plotting the limit points.  My blogging software ruined his code when he placed it in the comments section so I reproduce it here (but he’s also uploaded it to sagenb)

var('x malthus')
step(x,malthus) = malthus * x * (1-x)
stepfast = fast_callable(step, vars=[x, malthus], domain=RDF)

def logistic(m, step):
   # filter cycles
   v = .5
   for i in range(100):
       v = stepfast(v,m)
   points = []
   for i in range(100):
       v = stepfast(v,m)
       points.append((m+step*random(),v))
   return points

points=[]
step = 0.005
for m in sxrange(2.5,4,step):
   points += logistic(m, step)
point(points,pointsize=1).show(dpi=150)

One of the problems with the Linux MATLAB installer is that it doesn’t add any items to the menu in windows managers such as GNOME or KDE so you have to do it manually.  Now, there are several ways to do this but one of the easiest is to use the GUI tools provided with Ubuntu (assuming you are using Ubuntu of course).  For the record I am using Ubuntu 9.10 (Karmic Koala) with the default GNOME windows manager along with  MATLAB 2009b but I imagine that these instructions will work for a few other configurations too.

First things first, you’ll want to download a nice, scalable MATLAB icon since the ones included with MATLAB itself are a bit crude to say the least.  I tried to create one directly from MATLAB using the logo command and the Scalable Vector Graphics (SVG) Export of Figures package on The File Exchange but the result wasn’t very good (see below).

Fortunately someone called malte has created a much nicer .svg file (using inkscape) that looks just like the MATLAB logo and made it available via his blog.

Once you have the svg file you can start creating the shortcut as follows:

• On the GNOME Desktop click on System->Preferences->Main Menu
• Once the Main Menu program has started choose where you want your MATLAB icon to go and click on New Item.  In the screenshot below I have started to put it in the Programming submenu.

• Give your shortcut a name and put matlab -desktop in the command field. If you didn’t create a shortcut to the matlab executable when you installed the program then you may need to put the full path in instead (i.e. something like /opt/MATLAB/2009b/bin/matlab )

• Now let’s sort out the icon. Click on the spring on the left hand side of the above window to get the window below.

• Click on Browse and then select the folder that contains the .svg file you downloaded from malte’s blog.  I put it in my Desktop folder for the screenshot below.  Note that the svg file doesn’t appear in the preview pane.  Click on Open.

• Choose your icon.

• Click on Close and you are done

If you are using the Student version of MATLAB then the process is slightly different and has been covered by my friend Paul Brabban over at CrossedStreams.com.

Earlier today I was chatting to a lecturer over coffee about various mathematical packages that he might use for an upcoming Masters course  (note – offer me food or drink and I’m happy to talk about pretty much anything). He was mainly interested in Mathematica and so we spent most of our time discussing that but it is part of my job to make sure that he considers all of the alternatives – both commercial and open source. The course he was planning on running (which I’ll keep to myself out of respect for his confidentiality) was definitely a good fit for Mathematica but I felt that SAGE might suite him nicely as well.

“Does it have nice, interactive functionality like Mathematica’s Manipulate function?” he asked

Oh yes! Here is a toy example that I coded up in about the same amount of time that it took to write the introductory paragraph above (but hopefully it has no mistakes). With just a bit of effort pretty much anyone can make fully interactive mathematical demonstrations using completely free software. For more examples of SAGE’s interactive functionality check out their wiki.

Here’s the code:

def ftermSquare(n):
 return(1/n*sin(n*x*pi/3))

def ftermSawtooth(n):
 return(1/n*sin(n*x*pi/3))

def ftermParabola(n):
 return((-1)^n/n^2 * cos(n*x))

def fseriesSquare(n):
 return(4/pi*sum(ftermSquare(i) for i in range (1,2*n,2)))

def fseriesSawtooth(n):
 return(1/2-1/pi*sum(ftermSawtooth(i) for i in range (1,n)))

def fseriesParabola(n):
 return(pi^2/3 + 4*sum(ftermParabola(i) for i in range(1,n)))

@interact
def plotFourier(n=slider(1, 30,1,10,'Number of terms')
,plotpoints=('Value of plot_points',[100,500,1000]),Function=['Saw Tooth','Square Wave','Periodic Parabola']):
    if Function=='Saw Tooth':
     show(plot(fseriesSawtooth(n),x,-6,6,plot_points=plotpoints))
    if Function=='Square Wave':
     show(plot(fseriesSquare(n),x,-6,6,plot_points=plotpoints))
    if Function=='Periodic Parabola':
     show(plot(fseriesParabola(n),x,-6,6,plot_points=plotpoints))

I’ve just installed Ubuntu version 9.10 (Karmic Koala) from scratch on my new work machine and noticed that you can no longer right click on a file to encrypt it. The functionality is still there – it just isn’t available by default. To get it just install the seahorse-plugins package as follows

sudo apt-get install seahorse-plugins

I did a reboot to get the changes to take but there is probably a less drastic way of getting the job done.

I was recently given the task of converting a small piece of code written in R, the free open-source programming language heavily used by statisticians, into MATLAB which was an interesting exercise since I had never coded a single line of R in my life!  Fortunately for me, the code was rather simple and I didn’t have too much trouble with it but other people may not be so lucky since both MATLAB and R can be rather complicated to say the least.  Wouldn’t it be nice if there was a sort of Rosetta Stone that helped you to translate between the two systems?

Happily, it turns out that there is in the form of The MATLAB / R Reference by David Hiebeler which gives both the R and MATLAB commands for hundreds of common (and some not so common) operations.

While flicking through the 47 page document I noted that there are a few MATLAB commands for which David hasn’t found an R equivalent (possibly because there simply isn’t one of course).  For example, at number 161 of David’s document he describes the MATLAB command

yy=spline(x,y,xx)

which he describes as

‘Fit cubic spline with “not-a-knot” conditions (the first two piecewise cubics coincide,as do the last two), to points (xi , yi ) whose coordinates are in vectors x and y; evaluate at points whose x coordinates are in vector xx, storing corresponding y’s in yy.’

At the moment David doesn’t know of an R equivalent so if you are a R master then maybe you could help out with this extremely useful document?