## Archive for the ‘python’ Category

xkcd is a popular webcomic that sometimes includes hand drawn graphs in a distinctive style. Here’s a typical example

In a recent Mathematica StackExchange question, someone asked how such graphs could be automatically produced in Mathematica and code was quickly whipped up by the community. Since then, various individuals and communities have developed code to do the same thing in a range of languages. Here’s the list of examples I’ve found so far

- xkcd style graphs in Mathematica. There is also a Wolfram blog post on this subject.
- xkcd style graphs in R. A follow up blog post at vistat.
- xkcd style graphs in LaTeX
- xkcd style graphs in Python using matplotlib
- xkcd style graphs in MATLAB. There is now some code on the File Exchange that does this with your own plots.
- xkcd style graphs in javascript using D3
- xkcd style graphs in Euler

Any I’ve missed?

If you have an interest in mathematics, you’ve almost certainly stumbled across The Wolfram Demonstrations Project at some time or other. Based upon Wolfram Research’s proprietary Mathematica software and containing over 8000 interactive demonstrations, The Wolfram Demonstrations Project is a fantastic resource for anyone interested in mathematics and related sciences; and now it has some competition.

Sage is a free, open source alternative to software such as Mathematica and, thanks to its interact function, it is fully capable of producing advanced, interactive mathematical demonstrations with just a few lines of code. The Sage language is based on Python and is incredibly easy to learn.

The Sage Interactive Database has been launched to showcase this functionality and its looking great. There’s currently only 31 demonstrations available but, since anyone can sign up and contribute, I expect this number to increase rapidly. For example, I took the simple applet I created back in 2009 and had it up on the database in less than 10 minutes! Unlike the Wolfram Demonstrations Project, you don’t need to purchase an expensive piece of software before you can start writing Sage Interactions….Sage is free to everyone.

Not everything is perfect, however. For example, there is no native Windows version of Sage. Windows users have to make use of a Virtualbox virtual machine which puts off many people from trying this great piece of software. Furthermore, the interactive ‘applets’ produced from Sage’s interact function are not as smooth running as those produced by Mathematica’s Manipulate function. Finally, Sage’s interact doesn’t have as many control options as Mathematica’s Manipulate (There’s no Locator control for example and my bounty still stands).

The Sage Interactive Database is a great new project and I encourage all of you to head over there, take a look around and maybe contribute something.

My attention was recently drawn to a Google+ post by JerWei Zhang where he evaluates 2^3^4 in various packages and notes that they don’t always agree. For example, in MATLAB 2010a we have 2^3^4 = 4096 which is equivalent to putting (2^3)^4 whereas Mathematica 8 gives 2^3^4 = 2417851639229258349412352 which is the same as putting 2^(3^4). JerWei’s post gives many more examples including Excel, Python and Google and the result is always one of these two (although to varying degrees of precision).

What surprised me was the fact that they disagreed at all since I thought that the operator precendence rules were an agreed standard across all software packages. In this case I’d always use brackets since _I_ am not sure what the correct interpretation of 2^3^4 should be but I would have taken it for granted that there is a standard somewhere and that all of the big hitters in the numerical world would adhere to it.

Looks like I was wrong!

Late last year I was asked to give a talk at my University to a group of academics studying financial mathematics (Black-Scholes, monte-carlo simulations, time series analysis – that sort of thing). The talk was about the NAG Numerical Library, what it could do, what environments you can use it from, licensing and how they could get their hands on it via our site license.

I also spent a few minutes talking about Python including why I thought it was a good language for numerical computing and how to interface it with the NAG C library using my tutorials and PyNAG code. I finished off with a simple graphical demonstration that minimised the Rosenbrock function using Python, Matplotlib and the NAG C-library running on Ubuntu Linux.

“So!”, I asked, “Any questions?”

The very first reply was “Does your Python-NAG interface work on Windows machines?” to which the answer at the time was “No!” I took the opportunity to ask the audience how many of them did their numerical computing in Windows (Most of the room of around 50+ people), how many people did it using Mac OS X (A small group at the back), and how many people did it in Linux (about 3).

So, if I wanted the majority of that audience to use PyNAG then I had to get it working on Windows. Fortunately, thanks to the portability of Python and the consistency of the NAG library across platforms, this proved to be rather easy to do and the result is now available for download on the main PyNAG webpage.

Let’s look at the main differences between the Linux and Windows versions

**Loading the NAG Library**

In the examples given in my Linux NAG-Python tutorials, the cllux08dgl NAG C-library was loaded using the line

libnag = cdll.LoadLibrary("/opt/NAG/cllux08dgl/lib/libnagc_nag.so.8")

In Windows the call is slightly different. Here it is for CLDLL084ZL (Mark 8 of the C library for Windows)

C_libnag = ctypes.windll.CLDLL084Z_mkl

If you are modifying any of the codes shown in my NAG-Python tutorials then you’ll need to change this line yourself. If you are using PyNAG, however, then it will already be done for you.

**Callback functions**

For my full introduction to NAG callback functions on Linux click here. The main difference between using callbacks on Windows compared to Linux is that on Linux we have

C_FUNC = CFUNCTYPE(c_double,c_double )

but on Windows we have

C_FUNC = WINFUNCTYPE(c_double,c_double )

There are several callback examples in the PyNAG distribution.

That’s pretty much it I think. Working through all of the PyNAG examples and making sure that they ran on Windows uncovered one or two bugs in my codes that didn’t affect Linux for one reason or another so it was a useful exercise all in all.

So, now you head over to the main PyNAG page and download the Windows version of my Python/NAG interface which includes a set of example codes. I also took the opportunity to throw in a couple of extra features and so upgraded PyNAG to version 0.16, check out the readme for more details. Thanks to several employees at NAG for all of their help with this including Matt, John, David, Marcin and Sorin.

Over at Mathrecreation, Dan Mackinnon has been plotting two dimensional polygonal numbers on quadratic number spirals using a piece of software called Fathom. Today’s train-journey home project for me was to reproduce this work using Mathematica and Python.

I’ll leave the mathematics to Mathrecreation and Numberspiral.com and just get straight to the code. Here’s the Mathematica version:

- Download the Mathematica Notebook
- Download a fully interactive version that will run on the free Mathematica Player

p[k_, n_] := ((k - 2) n (n + 1))/2 - (k - 3) n; Manipulate[ nums = Table[p[k, n], {n, 1, numpoints}]; If[showpoints, ListPolarPlot[Table[{2*Pi*Sqrt[n], Sqrt[n]}, {n, nums}] , Joined -> joinpoints, Axes -> False, PlotMarkers -> {Automatic, Small}, ImageSize -> {400, 400}], ListPolarPlot[Table[{2*Pi*Sqrt[n], Sqrt[n]}, {n, nums}] , Joined -> joinpoints, Axes -> False, ImageSize -> {400, 400} ] ] , {k, 2, 50, Appearance -> "Labeled"} , {{numpoints, 100, "Number of points"}, {100, 250, 500}} , {{showpoints, True, "Show points"}, {True, False}} , {{joinpoints, True, "Join points"}, {True, False}} , SaveDefinitions :> True ]

A set of 4 screenshots is shown below (click the image for a higher resolution version). The formula for polygonal numbers, p(kn,n), is only really valid for integer k and so I made a mistake when I allowed k to be any real value in the code above. However, I liked some of the resulting spirals so much that I decided to leave the mistake in.

When Dan showed the images he produced from Fathom, he mentioned that he would one day like to implement these systems in a more open platform. Well, you can’t get much more open than Python so here’s a piece of code that makes use of Python and the matplotlib package (download the file – polygonal_spiral.py).

#!/usr/bin/env python import matplotlib.pyplot as plt from math import * def p(k,n): return(((k-2)*n*(n+1))/2 -(k-3)*n) k=13 polygonal_nums = [p(k,n) for n in range(100)] theta = [2*pi*sqrt(n) for n in polygonal_nums] r = [sqrt(n) for n in polygonal_nums] myplot = plt.polar(theta,r,'.-') plt.gca().axis('off') plt.show() |

Now, something that I haven’t figured out yet is why the above python code works just fine for 100 points but if you increase it to 101 then you just get a single point. Does anyone have any ideas?

Unlike my Mathematica version, the above Python program does not give you an interactive GUI but I’ll leave that as an exercise for the reader.

**If you enjoyed this post then you might also like:**

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!

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))

Just a quick note to say that I’ve released version 0.15 of PyNAG – my Python interface to the NAG C library and given it its own dedicated page on Walking Randomly. New features include

- Basic support for 165 functions from Mark 8 of the NAG C library.
- 29 simple but fully functional examples.
- Added initial support for NAG’s Complex structure.

If you have learned how to program then the odds are that you have written a program that displays a list of prime numbers. A prime number generator is a programming rite of passage that is almost up there with the ubiquitous ‘Hello World‘ application (for programmers who are going on to do numerical work at least).

There are various ways you might choose to do it but I bet that your first (or even your second or third) thought wouldn’t be ‘*I know – I’ll use a regular expression‘.* Someone known only as ‘Abigail’ obviously thought that the world desperately needed such a regex since he / she came up with the following incantation

^1?$|^(11+?)\1+$

Which is a regular expression that I assure you can be used to generate a list of all the prime numbers. If you are in the mood for a puzzle then sit back, grab a cup of coffee and try to work out HOW it can generate the primes.

If you give up then head over to Python Programming for a full explanation along with Python code that demonstrates that it really does work.

Can anyone come up with regular expressions for other interesting series of numbers?

**Update (27th August 2009): **Someone pointed me to an alternative explanation here.

I’ll be attending the 2009 EuroSciPy scientific Python conference in Leipzig next month – July 25th and 26th to be precise – and would love to hear from anyone else who’ll be attending.