A colleague recently sent me the following code snippet in R
> a=c(1,2,3,40) > b=a[1:10] > b  1 2 3 40 NA NA NA NA NA NA
The fact that R didn’t issue a warning upset him since exceeding array bounds, as we did when we created b, is usually a programming error.
I’m less concerned and simply file the above away in an area of my memory entitled ‘Odd things to remember about R’ — I find that most programming languages have things that look odd when you encounter them for the first time. With that said, I am curious as to why the designers of R thought that the above behaviour was a good idea.
Does anyone have any insights here?
Octave is a free, open source language for numerical computing that is mostly compatible with MATLAB. For years, the official Octave project has been strictly command-line only which puts many users off — particularly those who were used to the graphical user interface (GUI) of MATLAB. Unofficial GUIs have come and gone over the years but none of them were fully satisfactory for one reason or another.
As of the 3.8 release of Octave on 31st December 2013, this all changed when Octave started shipping with its own, official GUI. It is currently considered as ‘experimental’ by the Octave team and is obviously rough around the edges here and there but it is already quite usable.
The system includes
- An editor with syntax highlighting and code folding
- A debugger
- A file browser and workspace viewer
- The ability to hide and move different elements of the GUI around (e.g. you could swap the positions of the workspace and File Browser, or tear-off the editor into its own Window)
- A documentation browser
- The Octave command window
I’ve spent an hour or so playing with it today and like it a lot!
Thanks to Júlio Hoffimann Mendes who let me know about this new release.
I work for The University of Manchester where, among other things, I assist in the support of various high performance and high throughput computing systems. Exchanges such as the following are, sadly, becoming all too commonplace
Researcher: “Hi, I have an embarrassingly parallel research problem that needs a lot of compute resource. Can you help?”
Support: “Almost certainly, you could have access to our 2500 core Condor pool or maybe our 2000 core HPC system or any number of smaller systems depending on the department you are in. Let’s meet to discuss your requirements in more detail”
Researcher: “Sounds great. I am using [insert expensive commercial package here], could we install that on your systems?”
Support: “Not unless you pay a HUGE amount of money because you’ll need dozens or maybe hundreds of licenses. The licenses will cost more than our machines! Could you use [insert open source equivalent here] instead?”
Researcher: “A member of your team suggested that about 2 years ago but [insert expensive commercial package here] is easier to use, looks pretty and a single license didn’t seem all that expensive. It’ll take me ages to convert to [insert open source equivalent here]. Instead of splitting the job up and spreading it around lots of machines, can’t I just run it on a faster machine?”
Support: “Sorry but parallelism is the only real game in town when it comes to making stuff faster these days. I’m afraid that you’ll have to convert to [insert open source equivalent here], open your chequebook or wait until 2076 for your results to complete on your workstation.”
The moral of the story is that if you want your compute to scale, you need to ensure that your licensing scales too.
Apple make a big deal out of the fact that their app stores for iPhone and iPad contain thousands upon thousands of apps (or applications for relative oldies such as myself). Some of them are free of change, many of them cost money but I got to wondering how many of them were open source.
Possibly the the best list of iOS open source software I have found is available at maniacdev.com which, at the time of writing, includes 42 different applications complete with iTunes links and the all important links to source code. Another useful resource is open.iphonedev.com which is a regularly updated directory of open source apps and libraries for iOS. There’s some great stuff available including Battle for Wesnorth, SCI-15C Scientific calculator and TuxRider (based on Tux Racer).
Free as in Speech but not always Free as in Beer
One of the things you’ll notice about iOS open source apps is that they often cost money and sometimes quite a lot which is in stark contrast to what you may be used to. For example, Battle for Wesnorth can be had for no money at all on platforms such as Linux and Windows but the iPad version costs $5.99 at the time of writing. The more serious, SCI-15C Scientific calculator costs $19.99 right now which is rather steep for any iPhone app let alone an open source one.
Charging money for open source software may upset some people but doing so is usually not against the terms and conditions of the underlying license. The Free Software Foundation (inventors of the GPL, one of the most popular forms of open source license) has the following to say on the matter (original source)
“Free software” is a matter of liberty, not price. To understand the concept, you should think of “free” as in “free speech,” not as in “free beer.
Personally, I am happy to pay a few dollars for the iPad version of an open-source app if the developer has done a good job of the port. What does surprise me, however, is that it seems like no one has taken the source-code of these apps, recompiled them and then released free-of-charge versions on the app store. This wouldn’t be against the license conditions of licenses such as the GPL so why hasn’t it been done? I wouldn’t do it because I feel that it would be unfair to the developer of the iOS version but I would be surprised if everyone felt this way.
There are many open source applications that I’d love to see ported to iPad. Here’s my top three wants:
- FreeCiv (An Android port is on the way!)
- Gnuplot (This was done for Windows Mobile ages ago – see my review of it here)
- Octave (The iPad is more powerful than the laptop I did my PhD on. Plenty for basic Octave use)
Over to you….What do you think of the state of open source software on iPhone and iPad? Which applications would you most like to see ported? What are your favourite open source apps?
Update: 9th March 2011. Apparently, many of the open source applications currently available on the App store today violate the terms of licenses such as the GPL. The Inquirer has more details.
Christmas isn’t all that far away so I thought that it was high time that I wrote my Christmas list for mathematical software developers and vendors. All I want for christmas is….
- A built in ternary plot function would be nice
- Ship workbench with the main product please
- An iPad version of Mathematica Player
- Merge the parallel computing toolbox with core MATLAB. Everyone uses multicore these days but only a few can feel the full benefit in MATLAB. The rest are essentially second class MATLAB citizens muddling by with a single core (most of the time)
- Make the mex interface thread safe so I can more easily write parallel mex files
- More CUDA accelerated functions please. I was initially excited by your CUDA package but then discovered that it only accelerated one function (Matrix Multiply). CUDA accelerated Random Number Generators would be nice along with fast Fourier transforms and a bit more linear algebra.
- Release Mathcad Prime.
- Mac and Linux versions of Mathcad. Maple,Mathematica and MATLAB have versions for all 3 platforms so why don’t you?
- Produce vector versions of functions like g01bk (poisson distribution function). They might not be needed in Fortran or C code but your MATLAB toolbox desperately needs them
- A Mac version of the MATLAB toolbox. I’ve got users practically begging for it :)
- A NAG version of the MATLAB gamfit command
- A just in time compiler. Yeah, I know, I don’t ask for much huh ;)
- A faster pdist function (statistics toolbox from Octave Forge). I discovered that the current one is rather slow recently
- A Locator control for the interact function. I still have a bounty outstanding for the person who implements this.
- A fully featured, native windows version. I know about the VM solution and it isn’t suitable for what I want to do (which is to deploy it on around 5000 University windows machines to introduce students to one of the best open source maths packages)
- An Android version please. Don’t make it free – you deserve some money for this awesome Mathcad alternative.
- The fact that you give the Windows version away for free is awesome but registration is a pain when you are dealing with mass deployment. I’d love to deploy this to my University’s Windows desktop image but the per-machine registration requirement makes it difficult. Most large developers who require registration usually come up with an alternative mechanism for enterprise-wide deployment. You ask schools with more than 5 machines to link back to you. I want tot put it on a few thousand machines and I would happily link back to you from several locations if you’ll help me with some sort of volume license. I’ll also give internal (and external if anyone is interested) seminars at Manchester on why I think Spacetime is useful for teaching mathematics. Finally, I’d encourage other UK University applications specialists to evaluate the software too.
- An Android version please.
How about you? What would you ask for Christmas from your favourite mathematical software developers?
MATLAB contains a function called pdist that calculates the ‘Pairwise distance between pairs of objects’. Typical usage is
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.
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?
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.
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:
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:
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
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.
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
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
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
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
This time maxima won’t throw an error because syntactically you haven’t done anything wrong, you merely told it to do
as opposed to what you really wanted which is
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
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
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
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..)
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:
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.B,A.B);
which is another way of doing the same old
where here I’ve made a 2×2 matrix A, a 2×1 matrix B, with A, A denoting rows 1 and 2 of A respectively and matrix multiplied the rows of A by B (using A.B, A.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!