Archive for the ‘walking randomly’ Category

August 22nd, 2012

Someone recently asked me if WalkingRandomly.com had a facebook page.  Since it wasn’t much effort to create one, I have now done so.  I have no idea if this will be of any use to anyone but a first stab at it is at http://www.facebook.com/Walkingrandomly 

Now I have to decide what to do with it. Does anyone have any thoughts on a blog such as this having its own facebook profile?  Is it a good idea?  Will anyone make use of it or is it just pointless?

August 20th, 2012

My first post on WalkingRandomly was on 20th August 2007 and, to be honest, I had no idea what I was going to write about or if I would keep at it; I certainly didn’t think that I would still be going strong 5 years later.  Although I love to write, there is one reason why I keep on pounding out the posts here…you!

To all commenters, twitter followers, emailers and anyone else who has contacted me over the last 5 years regarding this blog….THANK YOU!  Without you, I would have given up years ago because no one wants to write reams of stuff that doesn’t get read.

There are various things one could write about for a post like this.  Maybe I could post reader statistics for example or discuss what I’ve got out of blogging over the last few years.  Alternatively, I could go down the route of posting links to some of the posts I’m most proud of (or otherwise!).  I considered all these things but realised that what was most important to me was just to say thank you to all of you.

Thanks!

August 20th, 2011

This little corner of the web that I like to call home has been online for 4 years as of today.  I’d like to thank everyone who reads and comments on the stuff I put up here, I hope you find my wanderings interesting and/or useful.

I’d also like to thank all of the people and organisations I have collaborated with over the years while writing articles here. I’m not going to mention names out of fear of missing someone out but you know who you are.  The only exception would be Matthew Haworth of Reason Digital who talked me into starting this up during his University of Manchester farewell party.  He also provides me with web hosting (surviving three slashdottings so far), technical support and, on occasion, lunch!  Without Matt, there would be no Walking Randomly.

Back on WRs first birthday I gave out some stats (143 RSS subscribers and 250 visits per day) and so I’ll so the same today (1645 RSS subscribers and 800 unique visits per day)

Finally, here is a list of some of my favourite articles from over the last 4 years.  I hope you enjoy reading them as much as I enjoyed writing them

June 8th, 2011

I’ve been a user of Ubuntu Linux for years but the recent emphasis on their new Unity interface has put me off somewhat.  I tried to like it but failed.  So, I figured that it was time for a switch to a different distribution.

I asked around on Twitter and got suggestions such as Slackware, Debian and Linux Mint.  I’ve used both Slackware and Debian in the past but, while they might be fine for servers or workstations, I prefer something more shiny for my personal laptop.

I could also have stuck with Ubuntu and simply installed GNOME using synaptic but I like to use the desktop that is officially supported by the distribution.

So, I went with Linux Mint.  It isn’t going well so far!

I had no DVDs in the house so I downloaded the CD version, burned it to a blank CD and rebooted only to be rewarded with

Can not mount /dev/loop0 (/cdrom/casper/filesystem.squashfs) on //filesystem.squashfs

I checked the md5sum of the .iso file and it was fine. I burned to a different CD and tried again. Same error.

I was in no mood for a trawl of the forums so I simply figured that maybe something was wrong with the CD version of the distribution – at least as far as my machine was concerned. So, I started downloading the DVD version and treated my greyhound to a walk to the local computer shop to buy a stack of DVDs.

When I got back I checked the .md5 sum of the DVD image, burned it to disk and…got the same error. A trawl of the forums suggests that many people have seen this error but no reliable solution has been found.

Not good for me or Linux Mint but at least Desmond (below) got an extra walk!

Desmond the greyhound

Update 1 I created a bootable USB memory stick from the DVD .iso to elimiate any problems with my burning software/hardware. Still get the same error message. MD5 checksum of the .iso file is what it should be:

md5sum ./linuxmint-11-gnome-dvd-64bit.iso
773b6cdfe44b91bc44448fa7b34bffa8  ./linuxmint-11-gnome-dvd-64bit.iso

My machine is a Dell XPS M1330 which has been running Ubuntu for almost 3 years.

Update 2: Seems that this bug is not confined to Mint. Ubuntu users are reporting it too. No fix yet though
https://bugs.launchpad.net/ubuntu/+bug/636711

Update 3: There is DEFINITELY nothing wrong with the installation media.  Both USB memory stick and DVD versions boot on my wife’s (much newer)HP laptop with no problem.  So, the issue seems to be related to my particular hardware.  This is like the good old days of Linux where installation was actually difficult.  Good times!

Update 4: After much mucking around I finally gave up on a direct install of Mint 11.  The installer is simply broken for certain hardware configurations as far as I can tell.  Installed Mint 10 from the same pen drive that failed for Mint 11 without a hitch.

Update 5: As soon as the Mint 10 install completed, I did an apt-get dist-upgrade to try to get to Mint 11 that way. The Mint developers recommend against doing dist-upgrades but I don’t seem to have a choice since the Mint 11 installer won’t work on my machine. After a few minutes I get this error

dpkg: error processing python2.7-minimal (--configure):
 subprocess installed post-installation script returned error exit status 3
Errors were encountered while processing:
 python2.7-minimal

This is mentioned in this bug report.  I get over that (by following the instructions in #9 of the bug report) and later get this error

p: cannot stat `/usr/lib/pango/1.6.0/module-files.d/libpango1.0-0.modules': No such file or directory
cp: cannot stat `/usr/lib/pango/1.6.0/modules/pango-basic-fc.so': No such file or directory
E: /usr/share/initramfs-tools/hooks/plymouth failed with return 1.
update-initramfs: failed for /boot/initrd.img-2.6.35-22-generic
dpkg: error processing initramfs-tools (--configure):
 subprocess installed post-installation script returned error exit status 1
Errors were encountered while processing:
 initramfs-tools

I fixed this with

sudo ln -s x86_64-linux-gnu/pango /usr/lib/pango

Trying the apt-get dist-upgrade again leads to

The following packages have unmet dependencies:
 python-couchdb : Breaks: desktopcouch (< 1.0) but 0.6.9b-0ubuntu1 is to be installed
 python-desktopcouch-records : Conflicts: desktopcouch (< 1.0.7-0ubuntu2) but 0.6.9b-0ubuntu1 is to be installed

Which, thanks to this forum post, I get rid of by doing

sudo dpkg --configure -a
sudo apt-get remove python-desktopcouch-records desktopcouch evolution-couchdb python-desktopcouch

A few more packages get installed before it stops again with the error message

Unpacking replacement xserver-xorg-video-tseng ...
Processing triggers for man-db ...
Processing triggers for ureadahead ...
Errors were encountered while processing:
 /var/cache/apt/archives/xserver-xorg-core_2%3a1.10.1-1ubuntu1.1_amd64.deb
E: Sub-process /usr/bin/dpkg returned an error code (1)

I get past this by doing

sudo apt-get -f install

Then I try apt-get upgrade and apt-get dist-update again…possibly twice and I’m pretty much done it seems.

Update 6: On the train to work this morning I thought I’d boot into my shiny new Mint system. However I was faced with nothing but a blank screen.  I rebooted and removed quiet and splash from the grub options to allow me to see what was going on. The boot sequence was getting stuck on something like checking battery state. Up until now I had only been using Mint while connected to the Mains. Well, this was the final straw for me.  As soon as I got into work I shoved in a Ubuntu 11.04 live disk which installed in the time it took me to drink a cup of coffee. I’ve got GNOME running and am now happy. 

My Linux Mint adventure is over.

April 11th, 2011

I love mathematics and I also love gadgets so you’d think that I’d be overjoyed to learn that there are a couple of new graphical calculators on the block.  You’d be wrong!

Late last year, Casio released the Prizm colour graphical calculator.  It costs $130 and its spec is pitiful:

  • 216*384 pixel display with 65,536 colours
  • 16Mb memory
  • The CPU is a SuperH 3 running at 58Mhz (according to this site)

Casio Prizm.

More recently, Texas Instruments countered with its color offering, the TI-NSpire CX CAS.  This one costs $162 (source) and its specs are also a bit on the weak side but quite a bit higher than the Casio.

  • 320*240 pixels with 65,536 colours
  • 100Mb memory
  • CPU? I have no idea.  Can anyone help?

If you are into retro-computing then those specs might appeal to you but they leave me cold. They are slow with limited memory and the ‘high-resolution’ display is no such thing. For $100 dollars more than the NSpire CX CAS I could buy a netbook and fill it with cutting edge mathematical software such as Octave, Scilab, SAGE and so on.  I could also use it for web browsing,email and a thousand other things.

I (and many students) also have  mobile phones with hardware that leave these calculators in the dust.  Combined with software such as Spacetime or online services such as Wolfram Alpha, a mobile phone is infinitely more capable than these top of the line graphical calculators.

They also only ever seem to be used in schools and colleges.  I spend a lot of time working with engineers, scientists and mathematicians and I hardly ever see a calculator such as the Casio Prizm or TI NSpire on their desks.  They tend to have simple calculators for everyday use and will turn to a computer for anything more complicated such as plotting a graph or solving equations.

One argument I hear for using these calculators is ‘They are limited enough to use in exams.‘  Sounds sensible but then I get to thinking ‘Why are we teaching a generation of students to use crippled technology?‘ Why not go the whole hog and ban ALL technology in exams?  Alternatively, supply locked down computers for exams that limit the software used by students.  Surely we need experts in useful technology, not crippled technology?

So, I don’t get it.  Why do so many people advocate the use of these calculators?  They seem pointless!  Am I missing something?  Comments welcomed.

Update 1: I’ve been slashdotted!  Check out the slashdot article for more comments.

Update 2: My favourite web-comic, xkcd, covered this subject a while ago.

Other posts you may find useful / interesting

March 17th, 2011

Part of my job at the University of Manchester is to help support the use of High Throughput Computing (HTC) services.  I am part of a team that works within Manchester’s Faculty of Engineering and Physical Sciences (Physics,Chemistry,Maths,Engineering,Computer Science,Earth Sciences) but we don’t just work with ‘our’ schools; we also collaborate with teams from all over the University along with the central Research Support team.

For example, say you are a University of Manchester researcher and have written a monte-carlo simulation in MATLAB where a typical run takes 5 hours to complete. In order to get good results you might want to run this simulation 1000 times which is going to take your desktop machine quite a while; about 7 months in fact! You can either wait for 7 months or you can contact us for assistance. We’ll then do the following for you:

  • Attempt to optimise your MATLAB code. Recent speed-ups range from 10% to 20,000%
  • Give you access to our Condor pool (currently peaking at 1600 processor cores, more expected in the near future)
  • Assist you in modifying your code and writing script wrappers to make best use of those 1600 cores.
  • If appropriate, put you in touch with colleagues who support traditional HPC (High Performance Computing) supercomputers, GPUs (Graphics Processing Units) and similar technology.

Once we’ve done our bit, you throw your code at our Condor pool and get your results in less than an evening! Obviously we are not confined to MATLAB; we’ve also assisted users of Python, Mathematica, C, FORTRAN, Amber, Morphy and many more.  Our job is to help researchers do research more quickly!

For a mathy science geek with a penchant for mathematical software it really doesn’t get much better than this!  I get to play with high end hardware and the latest software, I get to learn new areas of science and mathematics and I get to work with world-leading experts in a multitude of fields.  More importantly, I get to make a real difference.

Yep, this is one part of my job that I really love.  Expect to see more articles on High Throughput Computing and technologies such as Condor on WalkingRandomly in the very near future.

December 5th, 2010

A couple of years ago I wrote an article called Christmas gifts for math geeks and it has proven to be quite popular so I decided to write a follow up.  As I started thinking about what I might include, however, I started to realise that I had produced a list for science geeks instead.  So, here it is – my recommendations for gifts for the scientist in your life.

Mathematica 8 Home Edition – This is the full version of Mathematica, possibly my favourite piece of mathematical computer software, at the extremely low price of 195 pounds + VAT.  I know what you are thinking ‘Over 200 quid is not an extremely low price.’ and I would tend to agree.  It is, however, very good value since a commercial license costs several thousand pounds and Mathematica is as good as MATLAB with a whole slew of toolboxes.  Mathematica is possibly the most feature complete piece of mathematical software available today and is infinitely better than any dedicated graphical calculator.
Mathematica 8 logo

Bigtrak – I don’t have a Bigtrak but I used to have one back in the 1980s.  Is the science geek in your life into computers and 30-40 years old?  If so then there is a distinct possibility that their first foray into the world of computer programming was with a Bigtrak back when they were 8 or so – I mean, this thing can even do loops!  This isn’t identical to the original but it is a very close facsimile and would be great for budding computer nerds or their misty eyed old dad.

200-in-1 electronic project lab.  Now this one brings back fond memories for me since it was given to me for my 10th birthday and is probably the reason I studied physics at A-Level since A-Level physics included the study of basic electronics.  I did well in A-Level physics and enjoyed it so I chose theoretical physics for my degree later moving on to a PhD so you could argue that this piece of kit changed my life!

I was overjoyed when I discovered that it was still being sold and was immensely pleased when I received it as a birthday present once again when I was 28.

The first thing you need to know about this wonderful piece of kit is that it requires no soldering; you wire up all of the components using bendy little springs – nothing could be more simple.  There is also no need to be able to read schematic diagrams (although this can be a great way to learn how to) since each spring is numbered so producing your own AM radio transmitter can be as simple as joining spring 1 to spring 10 to spring 53 and so on.

The practical upshot of all of this is that you can approach this thing at a variety of levels.  In the first instance you can just have fun building and playing with the various circuits which include things like a crystal set radio, a Morse code transmitter, a light detector, a sound detector and basic electronic games.  Once you’ve got that out of your system you can start to learn the basics of electronics if you wish.

I have since discovered 300 in 1 and even 500 in 1 electronic project labs which look great and all but this is the one that will forever be in my heart.

Wonders of the Solar System – I have always loved (although never practised) astronomy and avidly followed the adventures of Voyagers 1 and 2 when I was small.  Since then, modern space probes such as Cassini-Huygens, Galileo and Mars Odyssey have added more to our knowledge of our astronomical backyard  and we now know a tremendous amount about the solar system.  In this series, Brian Cox of the University of Manchester takes us on a grand-tour around the solar system.  The imagery is fantastic, Cox’s enthusiasm is infectious and the science is awesome.  Yep, I quite like this DVD :)

2011 ‘Lightning Calculation’ calendar – Ron Doerfler writes a blog called Dead Reckonings that specialises in the lost arts of the mathematical sciences.  Last year he designed a 2010 Graphical Computing calendar and made the designs available for free to allow you to print your own.  Centred around ancient computing devices called nomograms, the calendar was beautiful and after Ron very kindly sent me a copy, I encouraged him to make a version that he could sell.  Well, I guess he took my advice because Ron is back with a 2011 calendar with the theme of ‘Lightning Calculations’ and this time he is selling it from Lulu.com.

Since Ron is an all round nice guy, he also offers a high resolution pdf of the calendar to allow you to print it off yourself but personally I plan on showing my support by putting an order in with Lulu.com.  Nice work Ron!

2011 math calendar

October 30th, 2010

Some time ago now, Sam Shah of Continuous Everywhere but Differentiable Nowhere fame discussed the standard method of obtaining the square root of the imaginary unit, i, and in the ensuing discussion thread someone asked the question “What is i^i – that is what is i to the power i?”

Sam immediately came back with the answer e^(-pi/2) = 0.207879…. which is one of the answers but as pointed out by one of his readers, Adam Glesser, this is just one of the infinite number of potential answers that all have the form e^{-(2k+1) pi/2} where k is an integer. Sam’s answer is the principle value of i^i (incidentally this is the value returned by google calculator if you google i^i – It is also the value returned by Mathematica and MATLAB). Life gets a lot more complicated when you move to the complex plane but it also gets a lot more interesting too.

While on the train into work one morning I was thinking about Sam’s blog post and wondered what the principal value of i^i^i (i to the power i to the power i) was equal to. Mathematica quickly provided the answer:

N[I^I^I]
0.947159+0.320764 I

So i is imaginary, i^i is real and i^i^i is imaginary again. Would i^i^i^i be real I wondered – would be fun if it was. Let’s see:

N[I^I^I^I]
0.0500922+0.602117 I

gah – a conjecture bites the dust – although if I am being honest it wasn’t a very good one.  Still, since I have started making ‘power towers’ I may as well continue and see what I can see. Why am I calling them power towers? Well, the calculation above could be written as follows:

small power tower

As I add more and more powers, the left hand side of the equation will tower up the page….Power Towers.  We now have a sequence of the first four power towers of i:

i = i
i^i = 0.207879
i^i^i = 0.947159 + 0.32076 I
i^i^i^i = 0.0500922+0.602117 I

Sequences of power towers

“Will this sequence converge or diverge?”, I wondered.  I wasn’t in the mood to think about a rigorous mathematical proof, I just wanted to play so I turned back to Mathematica.  First things first, I needed to come up with a way of making an arbitrarily large power tower without having to do a lot of typing.  Mathematica’s Nest function came to the rescue and the following function allows you to create a power tower of any size for any number, not just i.

tower[base_, size_] := Nest[N[(base^#)] &, base, size]

Now I can find the first term of my series by doing

In[1]:= tower[I, 0]

Out[1]= I

Or the 5th term by doing

In[2]:= tower[I, 4] 

Out[2]= 0.387166 + 0.0305271 I

To investigate convergence I needed to create a table of these. Maybe the first 100 towers would do:

ColumnForm[
 Table[tower[I, n], {n, 1, 100}]
 ]

The last few values given by the command above are

0.438272+ 0.360595 I
0.438287+ 0.360583 I
0.438287+ 0.3606 I
0.438275+ 0.360591 I
0.438289+ 0.360588 I

Now this is interesting – As I increased the size of the power tower, the result seemed to be converging to around 0.438 + 0.361 i. Further investigation confirms that the sequence of power towers of i converges to 0.438283+ 0.360592 i.  If you were to ask me to guess what I thought would happen with large power towers like this then I would expect them to do one of three things – diverge to infinity, stay at 1 forever or quickly converge to 0 so this is unexpected behaviour (unexpected to me at least).

They converge, but how?

My next thought was ‘How does it converge to this value?  In other words, ‘What path through the complex plane does this sequence of power towers take?”  Time for a graph:

tower[base_, size_] := Nest[N[(base^#)] &, base, size];
complexSplit[x_] := {Re[x], Im[x]};
ListPlot[Map[complexSplit, Table[tower[I, n], {n, 0, 49, 1}]],
 PlotRange -> All]

Power Tower Spiral

Who would have thought you could get a spiral from power towers? Very nice! So the next question is ‘What would happen if I took a different complex number as my starting point?’ For example – would power towers of (0.5 + i) converge?’

The answer turns out to be yes – power towers of (0.5 + I) converge to 0.541199+ 0.40681 I but the resulting spiral looks rather different from the one above.

tower[base_, size_] := Nest[N[(base^#)] &, base, size];
complexSplit[x_] := {Re[x], Im[x]};
ListPlot[Map[complexSplit, Table[tower[0.5 + I, n], {n, 0, 49, 1}]],
 PlotRange -> All]

Another power tower spiralThe zoo of power tower spirals

So, taking power towers of two different complex numbers results in two qualitatively different ‘convergence spirals’. I wondered how many different spiral types I might find if I consider the entire complex plane?  I already have all of the machinery I need to perform such an investigation but investigation is much more fun if it is interactive.  Time for a Manipulate

complexSplit[x_] := {Re[x], Im[x]};

tower[base_, size_] := Nest[N[(base^#)] &, base, size];

generatePowerSpiral[p_, nmax_] :=
  Map[complexSplit, Table[tower[p, n], {n, 0, nmax-1, 1}]];

Manipulate[const = p[[1]] + p[[2]] I;
 ListPlot[generatePowerSpiral[const, n],
 PlotRange -> {{-2, 2}, {-2, 2}}, Axes -> ax,
 Epilog -> Inset[Framed[const], {-1.5, -1.5}]], {{n, 100,
 "Number of terms"}, 1, 200, 1,
 Appearance -> "Labeled"}, {{ax, True, "Show axis"}, {True,
 False}}, {{p, {0, 1.5}}, Locator}]

Power Tower Manipulate
After playing around with this Manipulate for a few seconds it became clear to me that there is quite a rich diversity of these convergence spirals. Here are a couple more

Power Tower Spirals

Some of them take a lot longer to converge than others and then there are those that don’t converge at all:

non convergent power

Optimising the code a little

Before I could investigate convergence any further, I had a problem to solve: Sometimes the Manipulate would completely freeze and a message eventually popped up saying “One or more dynamic objects are taking excessively long to finish evaluating……”  What was causing this I wondered?

Well, some values give overflow errors:

In[12]:= generatePowerSpiral[-1 + -0.5 I, 200]

General::ovfl: Overflow occurred in computation. >>
General::ovfl: Overflow occurred in computation. >>
General::ovfl: Overflow occurred in computation. >>
General::stop: Further output of General::ovfl will be suppressed during this calculation. >>

Could errors such as this be making my Manipulate unstable?  Let’s see how long it takes Mathematica to deal with the example above

AbsoluteTiming[ListPlot[generatePowerSpiral[-1  -0.5 I, 200]]]

On my machine, the above command typically takes around 0.08 seconds to complete compared to 0.04 seconds for a tower that converges nicely; it’s slower but not so slow that it should break Manipulate.  Still, let’s fix it anyway.

Look at the sequence of values that make up this problematic power tower

generatePowerSpiral[-0.8 + 0.1 I, 10]

{{-0.8, 0.1}, {-0.668442, -0.570216}, {-2.0495, -6.11826},
{2.47539*10^7,1.59867*10^8}, {2.068155430437682*10^-211800874,
-9.83350984373519*10^-211800875}, {Overflow[], 0}, {Indeterminate,
  Indeterminate}, {Indeterminate, Indeterminate}, {Indeterminate,
  Indeterminate}, {Indeterminate, Indeterminate}}

Everything is just fine until the term {Overflow[],0} is reached; after which we are just wasting time. Recall that the functions I am using to create these sequences are

complexSplit[x_] := {Re[x], Im[x]};
tower[base_, size_] := Nest[N[(base^#)] &, base, size];
generatePowerSpiral[p_, nmax_] :=
  Map[complexSplit, Table[tower[p, n], {n, 0, nmax-1, 1}]];

The first thing I need to do is break out of tower’s Nest function as soon as the result stops being a complex number and the NestWhile function allows me to do this.  So, I could redefine the tower function to be

tower[base_, size_] :=
 NestWhile[N[(base^#)] &, base, MatchQ[#, _Complex] &, 1, size]

However, I can do much better than that since my code so far is massively inefficient.  Say I already have the first n terms of a tower sequence; to get the (n+1)th term all I need to do is a single power operation but my code is starting from the beginning and doing n power operations instead.  So, to get the 5th term, for example, my code does this

I^I^I^I^I

instead of

(4th term)^I

The function I need to turn to is yet another variant of Nest – NestWhileList

fasttowerspiral[base_, size_] :=
 Quiet[Map[complexSplit,
 NestWhileList[N[(base^#)] &, base, MatchQ[#, _Complex] &, 1,
 size, -1]]];

The Quiet function is there to prevent Mathematica from warning me about the Overflow error.  I could probably do better than this and catch the Overflow error coming before it happens but since I’m only mucking around, I’ll leave that to an interested reader.  For now it’s enough for me to know that the code is much faster than before:

(*Original Function*)
AbsoluteTiming[generatePowerSpiral[I, 200];]

{0.036254, Null}
(*Improved Function*)
AbsoluteTiming[fasttowerspiral[I, 200];]

{0.001740, Null}

A factor of 20 will do nicely!

Making Mathematica faster by making it stupid

I’m still not done though. Even with these optimisations, it can take a massive amount of time to compute some of these power tower spirals. For example

spiral = fasttowerspiral[-0.77 - 0.11 I, 100];

takes 10 seconds on my machine which is thousands of times slower than most towers take to compute. What on earth is going on? Let’s look at the first few numbers to see if we can find any clues

In[34]:= spiral[[1 ;; 10]]

Out[34]= {{-0.77, -0.11}, {-0.605189, 0.62837}, {-0.66393,
  7.63862}, {1.05327*10^10,
  7.62636*10^8}, {1.716487392960862*10^-155829929,
  2.965988537183398*10^-155829929}, {1., \
-5.894184073663391*10^-155829929}, {-0.77, -0.11}, {-0.605189,
  0.62837}, {-0.66393, 7.63862}, {1.05327*10^10, 7.62636*10^8}}

The first pair that jumps out at me is {1.71648739296086210^-155829929, 2.96598853718339810^-155829929} which is so close to {0,0} that it’s not even funny! So close, in fact, that they are not even double precision numbers any more. Mathematica has realised that the calculation was going to underflow and so it caught it and returned the result in arbitrary precision.

Arbitrary precision calculations are MUCH slower than double precision ones and this is why this particular calculation takes so long. Mathematica is being very clever but its cleverness is costing me a great deal of time and not adding much to the calculation in this case. I reckon that I want Mathematica to be stupid this time and so I’ll turn off its underflow safety net.

SetSystemOptions["CatchMachineUnderflow" -> False]

Now our problematic calculation takes 0.000842 seconds rather than 10 seconds which is so much faster that it borders on the astonishing. The results seem just fine too!

When do the power towers converge?

We have seen that some towers converge while others do not. Let S be the set of complex numbers which lead to convergent power towers. What might S look like? To determine that I have to come up with a function that answers the question ‘For a given complex number z, does the infinite power tower converge?’ The following is a quick stab at such a function

convergeQ[base_, size_] :=
  If[Length[
     Quiet[NestWhileList[N[(base^#)] &, base, Abs[#1 - #2] > 0.01 &,
       2, size, -1]]] < size, 1, 0];

The tolerance I have chosen, 0.01, might be a little too large but these towers can take ages to converge and I’m more interested in speed than accuracy right now so 0.01 it is. convergeQ returns 1 when the tower seems to converge in at most size steps and 0 otherwise.:

In[3]:= convergeQ[I, 50]
Out[3]= 1

In[4]:= convergeQ[-1 + 2 I, 50]
Out[4]= 0

So, let’s apply this to a section of the complex plane.

towerFract[xmin_, xmax_, ymin_, ymax_, step_] :=
 ArrayPlot[
 Table[convergeQ[x + I y, 50], {y, ymin, ymax, step}, {x, xmin, xmax,step}]]
towerFract[-2, 2, -2, 2, 0.1]

Rough Power Tower Fractal

That looks like it might be interesting, possibly even fractal, behaviour but I need to increase the resolution and maybe widen the range to see what’s really going on. That’s going to take quite a bit of calculation time so I need to optimise some more.

Going Parallel

There is no point in having machines with two, four or more processor cores if you only ever use one and so it is time to see if we can get our other cores in on the act.

It turns out that this calculation is an example of a so-called embarrassingly parallel problem and so life is going to be particularly easy for us.  Basically, all we need to do is to give each core its own bit of the complex plane to work on, collect the results at the end and reap the increase in speed. Here’s the full parallel version of the power tower fractal code

(*Complete Parallel version of the power tower fractal code*)

convergeQ[base_, size_] :=
 If[Length[
 Quiet[NestWhileList[N[(base^#)] &, base, Abs[#1 - #2] > 0.01 &,
 2, size, -1]]] < size, 1, 0];

LaunchKernels[];
DistributeDefinitions[convergeQ];
ParallelEvaluate[SetSystemOptions["CatchMachineUnderflow" -> False]];

towerFractParallel[xmin_, xmax_, ymin_, ymax_, step_] :=
 ArrayPlot[
 ParallelTable[
 convergeQ[x + I y, 50], {y, ymin, ymax, step}, {x, xmin, xmax, step}
, Method -> "CoarsestGrained"]]

This code is pretty similar to the single processor version so let’s focus on the parallel modifications.  My convergeQ function is no different to the serial version so nothing new to talk about there.  So, the first new code is

LaunchKernels[];

This launches a set of parallel Mathematica kernels. The amount that actually get launched depends on the number of cores on your machine.  So, on my dual core laptop I get 2 and on my quad core desktop I get 4.

DistributeDefinitions[convergeQ];

All of those parallel kernels are completely clean in that they don’t know about my user defined convergeQ function. This line sends the definition of convergeQ to all of the freshly launched parallel kernels.

ParallelEvaluate[SetSystemOptions["CatchMachineUnderflow" -> False]];

Here we turn off Mathematica’s machine underflow safety net on all of our parallel kernels using the ParallelEvaluate function.

That’s all that is necessary to set up the parallel environment.  All that remains is to change Map to ParallelMap and to add the argument Method -> “CoarsestGrained” which basically says to Mathematica ‘Each sub-calculation will take a tiny amount of time to perform so you may as well send each core lots to do at once’ (click here for a blog post of mine where this is discussed further).

That’s all it took to take this embarrassingly parallel problem from a serial calculation to a parallel one.  Let’s see if it worked.  The test machine for what follows contains a T5800 Intel Core 2 Duo CPU running at 2Ghz on Ubuntu (if you want to repeat these timings then I suggest you read this blog post first or you may find the parallel version going slower than the serial one).  I’ve suppressed the output of the graphic since I only want to time calculation and not rendering time.

(*Serial version*)
In[3]= AbsoluteTiming[towerFract[-2, 2, -2, 2, 0.1];]
Out[3]= {0.672976, Null}

(*Parallel version*)
In[4]= AbsoluteTiming[towerFractParallel[-2, 2, -2, 2, 0.1];]
Out[4]= {0.532504, Null}

In[5]= speedup = 0.672976/0.532504
Out[5]= 1.2638

I was hoping for a bit more than a factor of 1.26 but that’s the way it goes with parallel programming sometimes. The speedup factor gets a bit higher if you increase the size of the problem though. Let’s increase the problem size by a factor of 100.

towerFractParallel[-2, 2, -2, 2, 0.01]

The above calculation took 41.99 seconds compared to 63.58 seconds for the serial version resulting in a speedup factor of around 1.5 (or about 34% depending on how you want to look at it).
Power Tower Fractal

Other optimisations

I guess if I were really serious about optimising this problem then I could take advantage of the symmetry along the x axis or maybe I could utilise the fact that if one point in a convergence spiral converges then it follows that they all do. Maybe there are more intelligent ways to test for convergence or maybe I’d get a big speed increase from programming in C or F#?  If anyone is interested in having a go at improving any of this and succeeds then let me know.

I’m not going to pursue any of these or any other optimisations, however, since the above exploration is what I achieved in a single train journey to work (The write-up took rather longer though). I didn’t know where I was going and I only worried about optimisation when I had to. At each step of the way the code was fast enough to ensure that I could interact with the problem at hand.

Being mostly ‘fast enough’ with minimal programming effort is one of the reasons I like playing with Mathematica when doing explorations such as this.

Treading where people have gone before

So, back to the power tower story. As I mentioned earlier, I did most of the above in a single train journey and I didn’t have access to the internet. I was quite excited that I had found a fractal from such a relatively simple system and very much felt like I had discovered something for myself. Would this lead to something that was publishable I wondered?

Sadly not!

It turns out that power towers have been thoroughly investigated and the act of forming a tower is called tetration. I learned that when a tower converges there is an analytical formula that gives what it will converge to:

h(z) = -\frac{W(-ln z)}{ln(z)

Where W is the Lambert W function (click here for a cool poster for this function).  I discovered that other people had already made Wolfram Demonstrations for power towers too

There is even a website called tetration.org that shows ‘my’ fractal in glorious technicolor.  Nothing new under the sun eh?

Parting shots

Well, I didn’t discover anything new but I had a bit of fun along the way. Here’s the final Manipulate I came up with

Manipulate[const = p[[1]] + p[[2]] I;
 If[hz,
  ListPlot[fasttowerspiral[const, n], PlotRange -> {{-2, 2}, {-2, 2}},
    Axes -> ax,
   Epilog -> {{PointSize[Large], Red,
      Point[complexSplit[N[h[const]]]]}, {Inset[
       Framed[N[h[const]]], {-1, -1.5}]}}]
  , ListPlot[fasttowerspiral[const, n],
   PlotRange -> {{-2, 2}, {-2, 2}}, Axes -> ax]
  ]
 , {{n, 100, "Number of terms"}, 1, 500, 1, Appearance -> "Labeled"}
 , {{ax, True, "Show axis"}, {True, False}}
 , {{hz, True, "Show h(z)"}, {True, False}}
 , {{p, {0, 1.5}}, Locator}
 , Initialization :> (
   SetSystemOptions["CatchMachineUnderflow" -> False];
   complexSplit[x_] := {Re[x], Im[x]};
   fasttowerspiral[base_, size_] :=
    Quiet[Map[complexSplit,
      NestWhileList[N[(base^#)] &, base, MatchQ[#, _Complex] &, 1,
       size, -1]]];
   h[z_] := -ProductLog[-Log[z]]/Log[z];
   )
 ]

Tetration Manipulate
and here’s a video of a zoom into the tetration fractal that I made using spare cycles on Manchester University’s condor pool.

If you liked this blog post then you may also enjoy:

October 27th, 2010

I was having a quick browse through the current Top-500 supercomputer list and it appears to me that UK Universities hardly feature at all.  We’ve got Hector at Edinburgh University in at number 26 but that’s the UK’s National facility so we’d expect it to be good.

Other than that, we have The University of Southampton at Number 83, The University of Bristol at Number 378 and that seems to be about it (please correct me if I’m wrong).  There are a few more UK-based installations in the top-500 but none of them are in Universities as far as I can tell.

I’m a little surprised by this since I thought that we (UK academia) would give a better account of ourselves.  Does it matter that we don’t have more boxes in the top-500 within UK Academia?  Is everyone doing research on smaller scale Beowulf clusters, condor pools, GP-GPUs and fully tricked-out 12 core desktops these days?  Is the need for big boxes receding?  Any insight anyone?

April 13th, 2010

A shameless copy of The Endeavour’s weekend miscellany

Programming

Food

Technology

Life