Nano flower

A while ago I blogged about NanoArt 2007 , which is an online competition to create nanoart. Well, it turns out that I had completely forgotten about this competition, so I have done nothing more than form a few vague ideas about what I planned to create as a submission. The closing date for submissions is 31 December 2007, which leaves very little free time to create a worthwhile submission, so I wondered how I could create something using the bits and pieces that I happen to have lying around already.

One of the 3 electron microscope pictures that are supplied for you to start from is this

which they fittingly call "Nano Flower", and you are allowed to submit up to 5 "artistic" images for the competition.

I think I will choose to submit a 5-frame animation generated using a variant of the Belousov-Zhabotinsky reaction simulation software that I described here, where I use the "Nano Flower" image as an external input to steer the the reaction dynamics. This means that I am simulating a modified version of the BZ reaction, where the surface on which the reaction occurs is weakly contaminated so that the BZ reaction dynamics vary across the surface.

Here is a 5-frame animation that I quickly created, which fortuitously turns out to be cyclic because the period of the BZ dynamics is 5 frames.

Well, that's one frame of the animation. I gave up trying to get the animation to upload successfully to blogger.com, so I have put it on my web site here.

I have called this work "Fizzix" for fairly obvious reasons. Now I will wait to see whether I have broken the submission rules for the NanoArt 2007 competition, which don't mention animation as being an acceptable artform, not even 5-frame animations.

Update (7 January 2008):

All of the NanoArt 2007 competition entries can be viewed here, and the album containing my single competition entry is here (with the animation hosted here). My competition entry is a lot less colourful than other peoples' entries, so it is unlikely to attract much attention. I have to hope that people click down to my animation to discover the true nature of my competition entry.

Molecular manufacturing

The Center for Responsible Nanotechnology has a nice overview of its current findings on molecular manufacturing here, which summarises the pros and cons of being able to manipulate matter in a controlled way at the molecular scale.

Molecular manufacturing is fundamentally different from the "heat and stir" approach to building things with molecules, because it aims to directly control the placement and interaction of individual molecules. This opens up lots of new possibilities for building things (beneficial or dangerous), and the world will not be the same when the transformation to molecular manufacturing has occurred.

This revolution in manufacturing will make the Industrial Revolution seem trivial in comparison, and it will happen within a small number of decades (CRN says 2 decades). This means that it is likely to have a big impact on the lives of most people alive today, which gives you an incentive to read all about it here.

Merry Yuletide

Here is a movie in which I have created an animation of a festive texture on the surface of a trefoil knot. Merry Yuletide!

How did I create this movie?

I have been playing around with the Belousov-Zhabotinsky reaction simulation that I described here.

One thing that you can do is to find limit cycles of the simulation, where the state of the array of cells returns to a state that it visited earlier in the simulation. The whole sequence of states between two such repeats (including one of the end points) is then a limit cycle of the BZ simulation. Such limit cycles must exist because the state space is finite in size, so it it inevitable that the simulation must eventually revisit states that it visited earlier, though starting from a random initial state the likelihood that this occurs in a given timescale decreases rapidly as the size of the array increases.

One example that I particularly like is a cycle of length 10 that I found on a toroidal 13 by 13 array (using the same parameter values and colour scheme as here), and I show 2 of these cycles in the movie below:

This cycle is unusual because it has a low symmetry, and because it is pretty to look at despite its short length.

Using this cyclic BZ solution on a torus it is relatively easy to create the movie at the start of this posting, by using it to texture the toroidal surface of a trefoil knot, and choosing a colour scheme that maximises its festive feel.

Nerd test

I happened upon Nerdiometer, which describes an on-line test for assessing how nerdy you are. You can find the test here.

I couldn't resist trying this one. I suspect that the results of the test are biassed by a selection effect where the nerdier you are the more likely you are to take the test in the first place.

Anyway, here is my result:

17% of people score higher than me, and 82% score lower. I was disappointed to score so highly, because as I filled in the multiple-choice questionnaire I could see that the extremely nerdy answers were off the scale compared to where I stood. Maybe that was a trick to let me make the choices that I did without being too embarrassed about them.

The questionnaire was extremely selective in the areas it covered, because it focussed on asking about nerdy things to the almost complete exclusion of asking you about anything else. So, if you spend only small fraction of your time being a nerd, you will register highly on the nerd-scale as defined by this test.

I wish there was a way to find out how each answer was weighted to produce the overall result. Oh no! Wishing for that must in itself be a nerdy thing to do! Aargh! This is renormalisation gone mad!

Update (12 December 2007):

The result I got when I tried version 2 of the test (see here) is

OK, so now I am a "cool nerd" which sounds pretty good to me, but I am annoyed that my Science/Math and Technology/Computer scores are deemed to be so low! Who are the people who designed this test anyway? I'm not playing any more, maybe...

Belousov-Zhabotinsky reaction

Whilst exploring Wikipedia I came across the Belousov-Zhabotinsky reaction, and the phenomenology of its reaction dynamics such as the development of propagating waves of reactants is described here. Although I previously knew about the BZ reaction I had never implemented a simulation of it. It turns out to be remarkably easy to implement in Mathematica, and the BZ reaction dynamics can readily be displayed as a very pretty movie.

Although the original purpose of this type of simulation was to model the BZ reaction dynamics, you could imagine using variants of this type of simulation for generating interesting dynamic art forms.

A slightly modified form of an algorithm to simulate the BZ reaction due to Professor A K Dewdney is listed here.

(i) Select an integer q in the range 2 through 255. Cells may be in any of the states 1 through q.
(ii) Select two integers k1 and k2 in the range 1 through 8 and an integer g in the range 0 through 100.
(iii) In the transition from one "step" to the next the state of each cell is changed once according to rules (iv) - (vii).
(iv) A cell in state q changes to state 1.
(v) A cell in state 1 changes to state a/k1 + b/k2 + 1 where a is the number of neighbors of the cell which are in states 2 through q-1 and b is the number of neighbors in state q.
(vi) A cell in any of states 2 through q-1 changes to S/(9 - c) + g, where S is the sum of the states of the cell and its neighbors and c is the number of neighbors in state 1.
(vii) If the application of rule (v) or rule (vi) would result in a cell having a state > q then the state of that cell becomes q.

After some experimentation to find the fastest (and clearest) implementation in Mathematica, I created a function for implementing a BZ update step with toroidal boundary conditions:

update[state_, q_, k1_, k2_, g_] :=
Module[
{ones, qs, kernel, total, totalones, totalqs},
ones = Map[If[#==1, 1, 0]&, state, {2}];
qs = Map[If[#==q, 1, 0]&, state, {2}];
kernel = Table[1, {3}, {3}];
total = ListCorrelate[kernel, state, {{2,2},{2,2}}];
totalones = ListCorrelate[kernel, ones, {{2,2},{2,2}}] - ones;
totalqs = ListCorrelate[kernel, qs, {{2,2},{2,2}}] - qs;
MapThread[If[#>q, q, #]&[Switch[#1, 1, Floor[#4/k1+#3/k2+1], q, 1, _, Floor[#5/(9-#2)+g]]]&, {state, totalones, totalqs, 8-totalones-totalqs, total}, 2]
];

I used this BZ update function to generate the above movie of the BZ reaction dynamics, using suitable parameter values that I copied from here.

statesequence =
With[
{n=64, q=200, k1=2, k2=3, g=70},
NestList[update[#, q, k1, k2, g]&, RandomInteger[{1,q}, {n,Floor[4/3n]}], 250]
];
movie = Map[ArrayPlot[#, ColorFunction->"Rainbow"]&, statesequence];
ListAnimate[movie]

The poor image quality in the movie is due to the high degree of image compression, rather than due to any limitations of the above algorithm. Also, each frame is separately scaled and mapped onto the colour lookup table, which is probably not the best way of creating the movie.

Enigmatic comments

This is a placeholder for comments on postings at Enigmatics.

Please remember to identify which Enigma problem you are commenting on.

Enigmatic movements

I have decided that my Mathematica solutions to the New Scientist Enigma problems are cluttering up this blog, so I will move them to a separate blog called Enigmatics which I will generate using the incredibly useful WorkLife FrameWork; believe me, using WLFW is like having a 24/7 personal assistant organising all of your Mathematica notebooks. This allows me to present my solutions in a form that is much closer to the original Mathematica, but it does not allow for comments on blog entries.

I have looked into the possibility of managing the comments myself, but the effort needed to counter the inevitable spam is too great. I will instead provide a link in each Enigmatics blog entry back to this blog, which will allow comments to be made using the Blogger commenting system.

NanoArt 2007

Nanoart is defined here as:

NanoArt is a new art discipline related to micro/nanosculptures created by artists/scientists through chemical/physical processes and/or natural micro/nanostructures that are visualized with powerful research tools like Scanning Electron Microscope and Atomic Force Microscope.

The process of creating nanoart is described here as:

I bring the small world in front of my audience through high resolution electron microscope scans of natural micro or nanostructures and sculptures I create at micro or nano scale by physical or/and chemical processing. I take further steps by mixing the realistic images of this structures with abstract colors, digitally painting and manipulating the monochromatic electron scans, and finally printing them with long-lasting inks on canvas or fine art paper (giclee prints). This way, the scientific images become artworks and could be showcased for a large audience to educate the public with creative images that are appealing and acceptable.

In a nutshell, the electron microscope provides the raw image, and the artist colours it in to produce an enhanced (i.e. artistic) result. This could range from a trivial colour tinting of the raw image, through to an artistic rendition that is based only very loosely on the raw image. Of course, I favour the more artistic style of nanoart, and I have some ideas on new ways of creating such artwork, but there isn't room in this tiny blog posting to tell you all about it!

Anyway, back to the title of this posting: NanoArt 2007. This is an online competition to create nanoart, and the submission deadline is 31 December 2007. For people without ready access to an electron microscope to create their own raw images there are 3 high resolution monochromatic electron scans provided here (nano-flower), here (micro & nano), and here (nano-crystals). I think it would be fun to enter this competition.

Dresden art gallery

Here is a nice photo of a very small part of the world famous Dresden Art Gallery:

© Ben Bartels and André Schulz
http://www.chessbase.com/images2/2004/dresden/dresden018.jpg

Here is a snapshot that I took of the same scene (with some more context) in Second Life:

The exact location of this duplicate of the Dresden Art Gallery in Second Life is http://slurl.com/secondlife/Dresden%20Gallery/77/121/26. Visitors to SL do not need to pay anything just to explore.

I won't show you inside the Second Life duplicate of the art gallery, so you will have to visit it yourself to enjoy it. You can have a look at the WIRED report here if you want more details. Set aside a spare half-day for this adventure in SL. You won't be disappointed.

Sand painting

I found this video of sand painting here at Backreaction:

I think it is a wonderful art form, and it has given me some ideas ...

Enigma 1469

Here is how Mathematica can be used to solve New Scientist Enigma number 1469 in the 17 November 2007 issue.

WARNING: See the comments for details. There is an error in the solution given below that I need to fix. In working on a fix I discovered a bug in the HamiltonianCycle function in Mathematica's Combinatorica package (this bug has been acknowledged by the author of the Combinatorica package) which destroyed any chance of my fix working fully correctly, although the correct solution to Enigma 1469 could be seen lurking in amongst the erroneous clutter. I will abandon my solution to Enigma 1469 because it has absorbed too much of my time already, and the discovery of the HamiltonianCycle bug can stand as a testimonial to the cleverness of Enigma 1469. I suppose that I ought to leave the rest of this posting intact otherwise the comments would be orphaned, and also there are some useful techniques illustrated in my erroneous solution.

Load the Combinatorica package.
Needs["Combinatorica`"];

Derive the formula for triangular numbers.
triangularformula=Sum[i,{i,n}]
1/2 n (1 + n)

Compute a list of all triangular numbers less than 1000.
triangularnumbers = Table[triangularformula, {n,44}]
{1, 3, 6, 10, 15, 21, 28, 36, 45, 55, 66, 78, 91, 105, 120, 136, 153, 171, 190, 210, 231, 253, 276, 300, 325, 351, 378, 406, 435, 465, 496, 528, 561, 595, 630, 666, 703, 741, 780, 820, 861, 903, 946, 990}

Convert this to the corresponding list of lists of digits.
triangulardigits = IntegerDigits[triangularnumbers]
{{1}, {3}, {6}, {1, 0}, {1, 5}, {2, 1}, {2, 8}, {3, 6}, {4, 5}, {5, 5}, {6, 6}, {7, 8}, {9, 1}, {1, 0, 5}, {1, 2, 0}, {1, 3, 6}, {1, 5, 3}, {1, 7, 1}, {1, 9, 0}, {2, 1, 0}, {2, 3, 1}, {2, 5, 3}, {2, 7, 6}, {3, 0, 0}, {3, 2, 5}, {3, 5, 1}, {3, 7, 8}, {4, 0, 6}, {4, 3, 5}, {4, 6, 5}, {4, 9, 6}, {5, 2, 8}, {5, 6, 1}, {5, 9, 5}, {6, 3, 0}, {6, 6, 6}, {7, 0, 3}, {7, 4, 1}, {7, 8, 0}, {8, 2, 0}, {8, 6, 1}, {9, 0, 3}, {9, 4, 6}, {9, 9, 0}}

Construct the adjacency matrix of allowed adjacencies between triangular numbers.
adjacencies = Outer[(If[Last[#1]==First[#2], 1, 0])&, triangulardigits, triangulardigits, 1];

Convert the adjacency matrix to a set of rules that allow me to add a "tooltip" to each node of the corresponding the graph so that it can be explored by hovering the mouse over the graph: (1) extract the sparse array form of the adjacency matrix, (2) convert the sparse array into a list of rules for which graph nodes are linked by edges, (3) add a tooltip to each graph node.
adjacencies2 = SparseArray[adjacencies];
adjacencies3 = ((adjacencies2//ArrayRules//Most) /. HoldPattern[{i_,j_}->1] :> i->j);
adjacencies4 = adjacencies3 /. x_?IntegerQ :> Tooltip[x, triangularnumbers[[x]]];

Display the graph with tooltips.
GraphPlot[adjacencies4, DirectedEdges->True]
[Output omitted]

Convert the adjacency matrix to a graph. I need this representation in order to use graph manipulation algorithms.
g = FromAdjacencyMatrix[adjacencies, Type->Directed];

Extract the longest cycle in the graph, and extract the edges within this cycle.
longestcyclenodes = ExtractCycles[g]//Sort//Last//Most;
longestcycleedges = Partition[longestcyclenodes,2,1,{1,1}];

Display the graph with the longest cycle highlighted.
GraphPlot[adjacencies4, EdgeRenderingFunction -> (If[MemberQ[longestcycleedges,#2], {Red,Arrow[#1]}, {Opacity[0.3],LightGray,Arrow[#1]}]&)]

Convert the longest cycle to the corresponding a list of triangular numbers, and thence the sequence of digits encountered when going around the cycle.
longestcyclenumbers = triangularnumbers[[longestcyclenodes]];
longestcycledigits = longestcyclenumbers//IntegerDigits//Flatten
[Output omitted]

Compute the number of digits encountered when going around the longest cycle.
longestcycledigits//Length
[Solution omitted]

The Goldstone theorem

There is an interesting posting on Tommaso Dorigo's blog A Quantum Diaries Survivor entitled The Goldstone Theorem for Real Dummies. Quoting from his posting the Goldstone theorem is: "The spontaneous breaking of a continuous symmetry of the lagrangian generates massless scalars". Despite its claim to be for real dummies his derivation of the Goldstone theorem assumes familiarity with a lot of theory, which ensures that the "dumb" reader will not develop an intuitive feel for the masslessness of the particles (the Goldstone bosons, that is).

I am personally at ease with the level of theory that he uses, but I am in a very small minority in this respect. Wouldn't it be nice to write an intuitive description that does not assume that the reader has prior knowledge of much theory, but where the writer ensures that the intuitive description faithfully conforms to the underlying theory that is known to him. The Goldstone theorem can be understood intuitively without having to follow the detailed algebraic steps of its proof.

Here I will attempt to "prove" the Goldstone theorem using only words (i.e. no pictures, and no maths). Pictures would greatly improve the presentation, and maths would very succinctly summarise and generalise what is going on, but it is interesting to see what you can do with words alone.

Mass is defined as the propensity of the vector-valued field value alone (i.e. not including any field derivatives, or equivalently the field values at neighbouring points) to contribute to a restoring force that pushes the field towards an equilibrium value of the field. If within a local neighbourhood of field space there are direction(s) that have zero restoring force then each of these directions corresponds to a zero mass excitation.

It suffices to visualise a potential in field space whose (negative) gradient defines the field restoring force, and to ask ourselves what sort of shapes this potential could have in field space. We need to find places in field space where this potential has zero gradient, and then enquire whether the gradient of the potential is actually zero over the whole of a local neighbourhood rather than just at a single point, so that displacement of the field value within such a neighbourhood would encounter zero restoring force, and would thus correspond to a zero mass excitation.

The visualisation problem reduces to the classification of the stationarities of a scalar potential function in field space, which is something that is very easy to visualise (e.g. minima, maxima, saddle points, etc). The particular case that is relevant to the Goldstone theorem is when the potential function is constrained to be symmetric, e.g. under the continuous group of rotations about the origin of the vector-valued field. What sort of stationarities are allowed under such symmetry constraints? For a rotation symmetry the potential must be constant on each origin-centred spherical shell in field space, but different spherical shells have independent values of the potential. The potential is therefore fully described by its variation as a function of distance from the origin in field space (i.e. a 1-dimensional potential).

To find an equilibrium point you need to find a stationary point of the 1-dimensional version of the potential, then constancy of the full version of the potential over each spherical shell gives us zero gradient of the potential over the entire spherical shell that intersects this stationary point, which thus gives as many zero-mass degrees of freedom of the vector-valued field as there are directions within the spherical shell (i.e. one less than the dimensionality of the field). This is the Goldstone theorem.

Rereading what I have written above I am now uncertain whether it is of much use to anyone other than me. It is just the internal chatter (minus the pictures) that goes on inside my head when I think "Goldstone theorem", so it isn't guaranteed to be as intuitive for other people as it is for me.

Enigma 1468

I have had to skip a couple of Enigmas because there hasn't been time to do them. I may come back to them later, but don't hold your breath.

Here is how Mathematica can be used to solve New Scientist Enigma number 1468 in the 10 November 2007 issue. This solution draws heavily on the graph manipulation facilities in the Combinatorica package. Most of the work involves setting up an appropriate graph to represent the problem, after which the solution to the problem drops out by simply calling the HamiltonianCycle function.

Load the Combinatorica package for doing operations on graphs.

Needs["Combinatorica`"];

Define a function for computing 2D indexing from 1D indexing on a 6 by 6 array.

indices[i_] := {Quotient[i+5,6], Mod[i,6,1]};

Define a function for computing if a pair of nodes is adjacent in a Manhattan metric. You could use GridGraph[6,6] instead, but I want to show how to use an explicit adjacency function here.

manhattan[{i1_,j1_}, {i2_,j2_}] := Abs[i1-i2]==1&&j1==j2 || Abs[j1-j2]==1&&i1==i2;

Define a function for testing if adjacency between a pair of nodes is excluded.

exception[{i1_,j1_}, {i2_,j2_}] := Apply[Or, Map[{i1,j1}==#[[1]] && {i2,j2}!=#[[2]]&, {{{2,3},{1,3}}, {{3,4},{4,4}}, {{5,3},{6,3}}, {{5,6},{5,5}}}]];

Use these adjacency functions to define an adjacency matrix.

adjacency = SparseArray[{{_, _}?((manhattan[#1,#2] && Not[exception[#1,#2]])&[indices[#[[1]]], indices[#[[2]]]]&)->1}, {36,36}];

Use this adjacency matrix to define an adjacency function.

edge[from_, to_] := adjacency[[from, to]]==1;

Define the layout (i.e. embedding) of the graph.

coordinates = Table[index->{{0,1},{-1,0}}.indices[index], {index,36}];

Display the graph.

GraphPlot[g=MakeGraph[Range[36],edge], DirectedEdges->True, VertexLabeling->True, VertexCoordinateRules->coordinates]

[Output omitted]

Compute all of the Hamiltonian cycles in the graph. There is only one, as expected.

cycles = HamiltonianCycle[g,All]

[Output omitted]

Extract the edges in the Hamiltonian cycle.

cycleedges = Partition[cycles[[1]],2,1];

Define a function for drawing only the edges that lie on the Hamiltonian cycle.

cycleedge[coords_, vertices_, label_] := If[MemberQ[cycleedges,vertices], {Red,Thickness[0.01],Arrow[coords]}, {}];

Display the Hamiltonian cycle.

GraphPlot[MakeGraph[Range[36],edge], VertexLabeling->True, VertexCoordinateRules->coordinates, EdgeRenderingFunction->cycleedge]

[Output omitted]

Rotate the Hamiltonian cycle so that it starts at the correct node.

cycles2 = RotateLeft[#, Position[#,11][[1,1]]-1]&[Most[cycles[[1]]]]

[Output omitted]

Compute which positions in the Hamiltonian cycle fall on the nodes labelled "ENIGMA".

Map[Position[cycles2,#][[1,1]]&, Range[19,24]]

[Solution omitted]

Notebook computer woes

This is a very sad and nerdy posting, but those who have spent hours or days tracking down computer faults will know the state of mind that it creates, and the consequent need to drink large amounts of beer and write drivel!

I bought myself a new notebook computer (HP Pavilion dv9340ea) three months ago, and would you believe that a fault has already developed? Luckily, the fault turns out to be fairly benign so I can continue to use the computer, but the fault continues my unbroken record of always managing to buy notebook computers that go wrong.

What is the fault? The graphics display developed a habit of locking up with no warning, and the computer would then reboot showing various regular patterns on the boot screen, and then go to a "blue screen of death" (BSOD) informing me that nvlddmkm.sys had thrown error 116. The computer would then reboot again offering me the choice of booting into a "safe mode" of Windows Vista, which I duly selected, and I was then allowed back into the machine. At least I could get at any files that I wanted to extract from the computer.

So, off I went onto the internet searching for information about nvlddmkm.sys and error 116. In turns out that you can find lots of complaints about nvlddmkm.sys, which I found perversely gratifying because it meant that I had company. Various wise individuals offered solutions to this problem, which mainly focussed on manually uninstalling the old and buggy graphics driver (i.e. nvlddmkm.sys) and installing the latest version of this driver. I did this but it didn't fix the problem. I repeated this process several times, steadily enlarging the scope of my manual deletion of old files before installing new files, but I completely failed to fix the problem.

Maybe my nvlddmkm.sys problem was not in fact the same as the one that could be fixed by updating the graphics driver. Apparently, error 116 is quite a general fault that could be thrown by almost any error condition in nvlddmkm.sys.

I then did a full virus scan but it found nothing unusual.

I then decided to pay more attention to how my computer behaved without any special graphics driver installed. It showed the same regular patterns on the boot screen that preceded the BSOD that I mentioned at the start of this posting, but the BSOD itself did not occur, and I could boot into normal (i.e. not "safe mode") Windows Vista, so I had all Windows services available to me, but the screen resolution was relatively low because there was no special graphics driver installed.

The only other sign that something was wrong was a strange coloured granularity that showed up in what should have been featureless regions of colour on the screen. On rebooting, this fault occasionally disappeared, so I tried various experiments to check whether it was a thermal problem that went away temporarily if I switched off for a while. No luck there. I could not find a systematic way of influencing the fault.

In desperation I tried tapping the computer case with my finger, moving around all over the case, both on top and underneath. I noticed that when I tapped in one particular place under the case the strange coloured granularity disappeared at exactly the moment that I tapped the case. Jackpot! Although the fault quickly returned, and after a few repetitions of this cycle I eventually couldn't get the fault to go away at all, I had now absolutely convinced myself that my computer had a hardware fault that was probably in the graphics subsystem.

A loose memory chip maybe? The only memory that you could fiddle with was system RAM not video RAM, so reseating the system RAM modules didn't fix the problem, nor did I really expect it to.

However, I managed to find some photos of the system board of my notebook computer, and I noticed that immediately under the area of the case where I did my few successful taps was ... drum roll ... the graphics subsystem. I think this fairly convincingly demonstrates that something has gone wrong with the graphics hardware in my computer.

Now, all I have to do is to convince whoever it is that actually fixes my computer that they don't need to reformat my hard disk to check whether doing so cures the fault. In fact, it would be really nice if they would simply bring out a new system board to swap with my faulty one. What's the betting that things will not go as smoothly as this?

Update: After multiple emails and phone calls to HP they have agreed to service my computer under warranty, but they insist on restoring my hard disk to its factory settings. Resistance is useless! Since I don't relish going through the whole process of reinstalling all of my software, especially since there are lots of fiddly little things you have to do to get some of my old Windows XP software to work correctly under Windows Vista, I have decided to buy some disk mirroring software (i.e. Paragon Drive Copy 8.5). This software spent 2 1/2 hours mirroring my notebook computer's hard disk, and then announced that it had successfully done the job. Do I believe this? Do I believe that if it has done this part successfully, then it will also do the reverse mirroring successfully when I get my computer back from HP? I don't know what to think, but at least I have tried, and the Paragon software didn't cost too much (i.e. £30).

Update: Finally, after further phone calls to HP, they have agreed on specific dates to collect & return my notebook PC for its servicing-under-warranty. My main criticism of this process is not so much the large number of phone calls that I had to make, but the fact that HP uses call centres staffed by people who speak fluent English but with an impossible-to-understand accent. I had to ask several times for things to be repeated, and only with a certain amount of ingenuity and prior knowledge on my part could I understand what was being said. On more than one occasion I had to abort the phone call because I couldn't understand the accent at all, and on one occasion they aborted the phone call because of my inability to understand them!

Update: Finally! My notebook computer is back from the HP repair centre. The repair notes state that it had exactly the problem that I had originally diagnosed (i.e. a video hardware fault), and that the repair was a straight replacement of the motherboard. I hope that this fault was merely an "infant mortality", and that now I am likely to have years of trouble-free operation.

Enigma 1465

Here is how Mathematica can be used to solve New Scientist Enigma number 1465 in the 20 October 2007 issue.

This Enigma problem can be solved by using a very simple analytic approach which I give at the end of this posting. But first I give a much more flexible approach to illustrate the use of one of the Mathematica packages.

Solution 1: Brute force approach using the Combinatorica package.

Load the Combinatorica package.
Needs["Combinatorica`"]

Define a large enough grid graph size that the solution is contained within it.
n=11;

Build a directed grid graph which allows movement only rightwards and upwards. Note that the starting node is at the bottom-left hand corner. This type of directed structure forces all paths from the bottom left node to each other node to also be shortest paths.
g = MakeGraph[Range[n^2], Function[{x,y}, If[y==x+1 && Mod[x,n]>0 y==x+n, True, False]]]
-Graph:<220,121,directed>-

Compute the shortest path length from the bottom left corner node to each of the nodes in the graph. Because of the directed structure of the graph all paths are also shortest paths.
pathlengths = AllPairsShortestPath[g][[1]];

Generate a list of distinct shortest path lengths.
pathlengths2 = pathlengths//Union;

For each distinct shortest path length generate a list containing the number of distinct shortest paths to each node that has this shortest path length.
pathnumbers = Map[NumberOfKPaths[g,1,#]&, pathlengths2];

Extract the two cases where the same number of paths occurs 4 times. Since the solution lies in the quadrant opposite the starting node it must have a relatively large number of paths, so it can be selected by inspection of this result.
pathnumbers2 = Transpose[Cases[pathnumbers//Flatten//Tally, {_,4}]][[1]]
{10, 3003}

Solution 2: Quick approach using the known number of paths in a directed grid graph.

Compute a table each of whose entries is the number of paths from the top-left node of a directed (rightwards or downwards) grid graph.
With[{n=11}, Table[Binomial[i+j,j], {i,0,n-1}, {j,0,n-1}] // TableForm]
{
{1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1},
{1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11},
{1, 3, 6, 10, 15, 21, 28, 36, 45, 55, 66},
{1, 4, 10, 20, 35, 56, 84, 120, 165, 220, 286},
{1, 5, 15, 35, 70, 126, 210, 330, 495, 715, 1001},
{1, 6, 21, 56, 126, 252, 462, 792, 1287, 2002, 3003},
{1, 7, 28, 84, 210, 462, 924, 1716, 3003, 5005, 8008},
{1, 8, 36, 120, 330, 792, 1716, 3432, 6435, 11440, 19448},
{1, 9, 45, 165, 495, 1287, 3003, 6435, 12870, 24310, 43758},
{1, 10, 55, 220, 715, 2002, 5005, 11440, 24310, 48620, 92378},
{1, 11, 66, 286, 1001, 3003, 8008, 19448, 43758, 92378, 184756}
}

The three 6-fold paths mentioned in the first part of the problem can be seen in the correct positions in the table (i.e. (1,5), (5,1) and (2,2), where the top-left corner is (0,0)). The 4-times repeated solution(s) found above can be seen elsewhere in the table, and the larger solution lies in the lower right quadrant.

Ising model simulation - online software

I showed here how to use Mathematica (version 6) to create an interactive simulation of a 2-dimensional Ising model, but the software could be used only by people who owned their own copy of Mathematica. I reported here that Wolfram had now made available the Mathematica Player online conversion system for transforming Mathematica notebooks into a form that can be run by the Mathematica Player.

I have submitted my Ising model notebook for online conversion, and I have uploaded the converted notebook here (855 kB) plus a zipped version here (26 kB).

The Mathematica Player can be downloaded from here.

I would be interested to hear from anybody who does not have Mathematica (version 6) whether they can successfully get this converted Ising model notebook to work. I would particularly like comments on the ease of use of the zipped version, because this is by far the more disk-space-efficient way of doing things.

Update: Since no-one has responded yet to say that that they have succeeded in getting my converted Ising model notebook to work with Mathematica Player, I decided to install Mathematica Player on my old notebook PC (which didn't have Mathematica already installed) to check things out for myself. I downloaded the zipped version of the converted Ising model notebook, and fired it up in Mathematica Player. The contents of the converted notebook look perfectly OK, but seem to be completely inactive, so you can't run the Ising model simulation. I haven't got the faintest idea what is wrong here, so I'll now go away to try to sort this problem out.

Mathematica Player online conversion

At last! We can now use Mathematica to write software that we can publish on our websites and which can be run by anybody, i.e. this is not limited to end-users who own their own copy of Mathematica. This development was quietly announced in a Wolfram blog posting entitled The Day that Documents and Applications Merged; at least, that's where I learnt about it. It takes the form of a Mathematica Player online conversion system for transforming Mathematica notebooks into a form that can be run by the Mathematica Player. The blog posting summarises this as:

But the new online conversion service opens up almost any notebook. If you're a Mathematica user, that means you don't have to share your work just with other users anymore. Create your dynamic applications and share them with anyone you'd like. You can even post them on your website and link people to the Player download - a great way to publish research papers, books, or interesting computational models.

This changes everything! Long-time users of Mathematica such as me have been gnashing their teeth for 20 years at the fact that, although they have a uniquely powerful tool in Mathematica, any software that they write has been limited to being used by only those people who have their own copy of Mathematica. Unfortunately, there is almost zero chance that any randomly selected person owns Mathematica.

Many years ago I made the decision to use Mathematica because it was my "secret weapon" that allowed me to do things that would be much more cumbersome to do with other software tools such as C++, Matlab, etc. It allowed me to do a lot more research than would otherwise have been the case. In fact, I covered ground so fast that people couldn't keep up with me; that's praise for Mathematica rather than an inflated view of my ability to do research. I have received almost universal criticism from those around me (who use C++, Matlab, etc), because the fruits of my Mathematica labours are not available for them to use, at least not without a lot of effort to translate the software into a form that satisfies them (i.e. into C++, Matlab, etc). My use of Mathematica is regarded by these people as a selfish act, because it makes life easy for me and difficult for them. I can see their point! However, I am not one to be swayed by other peoples' misguided opinions, so I carried on regardless.

Of course, the Mathematica Player online conversion system comes with some restrictions, which limit its use to one of the following categories:

1. Educational purposes in a not-for-profit, educational institution.
2. Publicly accessible or published research content.
3. Demonstration purposes, including at commercial or governmental organizations.
4. None of the above. I am interested in commercial distribution.

This set of alternatives seems remarkably generous to me, after 20 years of not being able to use Mathematica in this way. The category above that particularly interests me is (2), which allows me to include working software with research papers that I publish. In fact, I wonder whether I will move to a model in which I seamlessly combine software and research paper in the Mathematica notebook that I feed to the Mathematica Player online conversion system (I assume this is possible).

Thank you so much Wolfram! I hope there aren't hidden caveats that break what appears to be a fully working model.

Random sculptures

Here is a piece of Mathematica code for generating interesting random sculptures by generating 2-dimensional manifolds embedded in 3-dimensional space.

Define a function that computes a random sum of sinewaves, where the maximum frequency is controlled by the n parameter. Differently seeded versions of this function will be used to determine the displacement of the 2-dimensional manifold in each of the directions in 3-dimensional space.

f[x_,n_]:=Sum[RandomReal[{-1,1}]Cos[i x+RandomReal[2\[Pi]]], {i,0,n}];

Display a random sculpture. The n parameter controls the complexity of the sculpture.

With[{n=1}, ParametricPlot3D[Table[f[u,n]f[v,n],{3}]//Evaluate, {u,0,2\[Pi]},{v,0,2\[Pi]}, PlotStyle->Opacity[0.25], Mesh->False, Boxed->False, Axes->False]]

Here are 3 examples of the sort of result that you can obtain. Not bad for a simple piece of code!

Enigmatics

I thought that I should explain the purpose of me posting Mathematica solutions to the New Scientist Enigma problems.

The main purpose is to show how one can quickly get a solution to each problem (usually by brute force attack) using the unique processing abilities of Mathematica. Each posting shows the first line of reasoning that led me to the solution, so my reasoning is usually not elegant and is not intended to be elegant; the goal is to find the correct solution and to do it quickly. Sometimes this means that I deliberately use an unusual sequence of manipulations in order to illustrate a particular style of programming. Occasionally this means that I use a peculiar Mathematica idiom that I just happen to know.

Actually, I do worry about the elegance of solutions, but optimising the elegance of a solution is completely different from simply finding the correct solution to a problem. In the case of the New Scientist Enigma problems I strongly suspect that most of them are constructed in a brute force way (i.e. start with a large set of potential solutions, then impose various perverse constraints, until there is only one member of the set remaining), so there is no guarantee that there are slick and elegant solutions to them.

Enigma 1464

Here is how Mathematica can be used to solve New Scientist Enigma number 1464 in the 13 October 2007 issue.

Generate a list of leap years in the 20th century during which the births can occur. The year 1900 is excluded because it was not a leap year.
lpyrs = Table[i, {i, 1904, 2000, 4}]
{1904, 1908, 1912, 1916, 1920, 1924, 1928, 1932, 1936, 1940, 1944, 1948, 1952, 1956, 1960, 1964, 1968, 1972, 1976, 1980, 1984, 1988, 1992, 1996, 2000}

Generate a list of all days during the above leap years in {yyyy, mm, dd} format.
(days = Flatten[Map[Table[Take[DateList[t],3], Evaluate[{t,#,#+365*86400,86400}&[AbsoluteTime[{#}]]]]&, lpyrs], 1]) // Length
9150

Convert this list of days into DDMMYY format.
(days2 = Map[(10000#[[3]] + 100#[[2]] + If[#[[1]]>=2000, #[[1]]-2000, #[[1]]-1900])&, days]) // Length
9150

Select the cases which are square days.
sqrdays=Select[days2, IntegerQ[Sqrt[#]]&]
{300304, 10404, 200704, 40804, 121104, 91204, 60516, 10816, 110224, 101124, 20736, 11236, 150544, 131044, 80656, 70756, 20164, 50176, 30276, 30976, 51076, 240100, 260100, 230400, 270400, 220900, 280900}

Split this list into sublists for each year.
sqrdays2 = Split[sqrdays, Take[IntegerDigits[#1],-2]==Take[IntegerDigits[#2],-2]&]
{{300304, 10404, 200704, 40804, 121104, 91204}, {60516, 10816}, {110224, 101124}, {20736, 11236}, {150544, 131044}, {80656, 70756}, {20164}, {50176, 30276, 30976, 51076}, {240100, 260100, 230400, 270400, 220900, 280900}}

Select the cases that have 6 or more different square days.
sqrdays3 = Select[sqrdays2, Length[#]>=6&]
{{300304, 10404, 200704, 40804, 121104, 91204}, {240100, 260100, 230400, 270400, 220900, 280900}}

Split the list of square days into sublists for each month.
sqrdays4 = Table[Select[sqrdays, FromDigits[Take[IntegerDigits[#],{-4,-3}]]==i&], {i,12}]
{{20164, 50176, 240100, 260100}, {110224, 30276}, {300304}, {10404, 230400, 270400}, {60516, 150544}, {80656}, {200704, 20736, 70756}, {40804, 10816}, {30976, 220900, 280900}, {131044, 51076}, {121104, 101124}, {91204, 11236}}

Select from all of the cases with 6 or more square days only those for which there is but a single square day in that month. These are the cases where the birth and death days are forced to lie in different months. There is one case that satisfies this criterion, which is thus the grandfather's date of birth.
Select[sqrdays3//Flatten, Length[sqrdays4[[Position[sqrdays4,#1][[1,1]]]]]==1&]
{300304}

Geodesic dome construction

I have been thinking about how to automate the design of geodesic domes using Mathematica. I started off by looking at what Google Earth could see of the geodesic domes at the Eden Project (see 50.3612°N, -4.74483°W in Google Earth). Here is a clear view of the structure of one of the domes.

There is a central pentagon at the highest point of the dome, which is surrounded by 5 hexagons each of which is divided into 6 small triangles, which are in turn surrounded by hexagons (which are not divided into triangles). If you look very carefully you will find other pentagons lower down the dome.

How could this geodesic dome be designed in a semi-automatic way? I have worked out a way of doing this using Mathematica which goes as follows.

Load the library of polyhedron operations so that you can derive tessellated versions of polyhedra.

Needs["PolyhedronOperations`"];

Fix the order of the tessellation. A larger order will eventually give you a geodesic dome with smaller faces.

n=6;

Display an icosahedron and its tessellation, and extract the faces of the tessellation for later use.

GraphicsRow[Map[Graphics3D,{#, faces=Geodesate[#,n]}&[PolyhedronData["Icosahedron", "Faces"]]]]

The main problem now is to find a way of extracting pentagons and hexagons from all of the triangles that appear in the above tessellation. This is not trivial because it has to be done in such a way that the extracted polygons neither overlap nor have any gaps between them. The strategy that I use to achieve this is to gradually build up a set of polygon centres in such a way that each new centre is exactly far enough away from all of the centres found thus far that it defines a new polygon that neither overlaps nor leaves a gap with the polygons found thus far.

Define a function for computing all of the triangle face edges surrounding a set of vertices.

faceedges[perimeters_, centres_] := Flatten[Map[Cases[perimeters, {x___,#,y___} :> {x, y}]&, centres], 1] // Union;

Define a function for computing all of the vertices surrounding a set of triangle face edges.

facecentres[perimeters_, edges_] := Map[Cases[perimeters, {x___,#[[1]],y___,#[[2]],z___} | {x___,#[[2]],y___,#[[1]],z___} :> {x,y,z}]&, edges] // Flatten // Union;

Iterating the above pair of functions causes the set of vertices to grow until it covers the whole of the original set with pentagon and hexagon centres. As long as the order of the tessellation is chosen appropriately, then this algorithm will give the centres needed to construct a geodesic dome. Have a play to see what the algorithm does for you, and how it might be improved.

Fix the starting vertex to be the one at the top of the tessellated icosahedron.

startingcentre = {Position[faces[[1]], {0,0,1}][[1,1]]};

Derive the pentagon and hexagon centres and edges of the geodesic dome.

centres = Nest[
(
faceedgesold = faceedges[faces[[2,1]], #];
centresnew = facecentres[faces[[2,1]], faceedgesold]
)&,
startingcentre,
5
];
edges = faceedges[faces[[2,1]], centres];

Display the geodesic dome.

Graphics3D[{
{Opacity[0.75], faces},
{PointSize[0.02], Red, Map[Point[faces[[1,#]]]&, centres]},
{Thickness[0.01], Blue, Map[Line[{faces[[1,#[[1]]]], faces[[1,#[[2]]]]}]&, edges]}
}
]

If you compare this result with the photo of the geodesic dome at the Eden Project then you can see that it has all of the essential features correct. There is a pentagon at the top of the dome, which is surrounded by hexagons, and lower down the dome there are more pentagons.

Enigma 1463

Here is how Mathematica can be used to solve New Scientist Enigma number 1463 in the 6 October 2007 issue. Because this problem is so easy to solve I don't bother to make use of symmetry to simplify the search for a solution.

Generate a list of permutations of the digits 1, 2, ..., 9.
(digitperms = Permutations[Range[9]]) // Length
362880

Arrange each of these sets of 9 digits into a 3 by 3 array and generate the corresponding 12 3-digit numbers by reading rows and columns both forwards and backwards.
(numbers = Map[Flatten[Map[{FromDigits[#], FromDigits[Reverse[#]]}&, Join[Partition[#,3], Transpose[Partition[#,3]]]], 1]&, digitperms]) // Length
362880

Select the cases which satisfy the required divisibility conditions.
(numbers2 = Select[numbers, Apply[And, Table[Mod[Count[Map[Mod[#,i]&, #], 0], i]==0, {i,9}]]&]) // Length
8

Extract the minimum and maximum 3-digit numbers from these (symmetrically related) results.
{Min[#], Max[#]}&[numbers2]
{164, 971}

Rotation group topology

Here is a video animation that illustrates the topology of the rotation group SO(3), which is the group of rotations in 3-dimensional space. My apologies for the small size and poor visual quality of the video - I am still experimenting with how best to generate animations and to upload them to blogger.

A rotation in 3-dimensional space can be done about an axis pointing in any direction, so there is a whole spherical surface of possible axis directions, which is represented by the sphere in the video. Each point on this spherical surface represents one possible orientation of the rotation axis, and the whole spherical surface is a 2-dimensional manifold representing all of the possible rotation axes. For the purpose of this blog posting only the topology of this manifold is important (i.e. it has a spherical topology); the detailed geometry of the surface is not important.

Once you have decided where the rotation axis points you then rotate about this axis, which generates a 1-dimensional angle subspace that lives outside the 2 angle dimensions that already live within the spherical manifold. Thus far the overall manifold is the original 2-dimensional spherical manifold (generated by the choice of axis direction) times a 1-dimensional manifold (generated by the choice of rotation). If you rotate by one full circuit (i.e. the rotation angle varies from 0 to 2Pi) then you appear to return to the starting point in the 1-dimensional angle subspace, so it seems that this subspace has a standard circular topology. But not all is what it seems to be.

The spherical manifold of possible rotation axes has a hidden degeneracy, because each rotation axis appears twice over. Each point on the sphere and its opposite point (i.e. its antipode) define rotation axes that are parallel, except that they have the opposite physical effect when a given rotation angle is applied about each of these two axes, i.e. the rotation axes are anti-parallel. So a rotation of x about one axis produces the same physical effect as a rotation of -x about its antipodal axis.

How does this degeneracy manifest itself topologically? The additional 1-dimensional angle subspace that is attached to each point on the spherical manifold must be identified with (i.e. glued to) the corresponding 1-dimensional angle subspace that is attached to the antipodal point on the spherical manifold.

This gluing-together of what were hitherto separate parts of the overall manifold imposes the constraint that a rotation of x about one axis produces the same physical effect to a rotation of -x about its antipodal axis. These constraints impose additional topological structure on the original 2-dimensional spherical manifold (generated by axis directions) times a 1-dimensional manifold (generated by rotations), because now antipodal points on the 2-dimensional spherical manifold are glued together via their attached 1-dimensional rotation angle subspaces.

This final glued version of the overall manifold captures all of the topological structure of the group of rotations in 3-dimensional space (i.e. SO(3)). It does not have the same topology as the unglued version because the gluing causes what were hitherto separate points on the manifold to become identified.

This gluing process is illustrated in the video, which shows how the gluing can be implemented by gradually morphing the 1-dimensional manifolds living at antipodal points on the 2-dimensional spherical manifold, until they eventually lie alongside each other ready to be glued together, and then they are finally coalesced into a single 1-dimensional manifold. Colour helps to identify corresponding points on the pair of yet-to-be-glued 1-dimensional manifolds, so that when the colours align with each other then the gluing operation can take place. It is clear from the video that this forces the pair of 1-dimensional manifolds to lie along a diameter of the sphere. Remember that the interior of the sphere is hitherto unused, so it is free to be used to represent rotation angle. Following the original 1-dimensional circular path from a point on the sphere back to itself again now corresponds to following a 1-dimensional path along a diameter of the sphere to its antipodal point, and then instantaneously returning to the starting point which has been glued to its antipodal point.

What consequences does this non-trivial topology have for rotations in 3-dimensional space? A rotation through 2Pi now corresponds to tracing a path along the diameter of the sphere from a point on the surface of the sphere to its antipodal point, and then returning instantaneously to the starting point. There is a non-trivial loop here because the path has clearly travelled somewhere else (i.e. the antipodal point) before returning home, and the loop can't be made to go away by moving it around inside (or on the surface of) the sphere, because the constraint that a pair of points on the path are antipodes is locked into the definition of the path.

So rotating by 2Pi in 3-dimensional space does not return you to where you started, because the path that you follow has a persistence that is enforced by the above antipodal locking-in effect. Of course, the physical object that you rotate has to have an appropriate physical sensitivity to rotation in order to observe this effect. A rotational scalar (e.g. the temperature of an object) will not show this effect, nor will a rotational vector (e.g. directions in 3-dimensional space) show this effect. An example of an object that does show this and other related effects is the Dirac Belt, where an object (i.e. the belt buckle) is rotated whilst the belt itself is used to retain a memory of what rotations have been applied to the buckle. After 2Pi rotation the belt is twisted so the overall state of the system is not the same as before the rotation was applied. The key to making an object that is physically sensitive to a rotation of 2Pi is that it must retain a memory of what rotations have been applied, which is equivalent to the object knowing about the path through the sphere discussed above. Obviously, there are many types of objects that exhibit this 2Pi sensitivity.

The best bit of all is that the locking-in effect of the pair of antipodal points can easily be neutralised by making an additional rotation of 2Pi, so that the overall rotation is 4Pi. This is illustrated in the case of the Dirac Belt, where the twist of the belt is zero after a 4Pi rotation, so the belt is not physically sensitive to rotating its buckle by 4Pi. The reason that the antipodal locking-in effect goes away when you rotate by 4Pi (i.e. make two circuits along a diameter through the sphere in the video) is that there are now two pairs of locked-in antipodal points, which can separately be moved around inside (or on the surface of) the sphere in such a way that the overall length of the loop can be made to go to zero, i.e. equivalent to no path at all. This interesting manipulation is not illustrated in the video.

So the topology of the rotation group SO(3) (i.e. the group of rotations in 3-dimensional space) is non-trivial as shown in the video, and it has highly non-trivial physical consequences such as those illustrated by the Dirac Belt.

Ig Nobel Prizes 2007

The Ig Nobel Prizes 2007 have now been awarded. Quoting from News24 here this year's winners include:

1. Chemistry: Mayu Yamamoto of the International Medical Centre of Japan, for developing a way to extract vanillin, or vanilla fragrance and flavouring, from cow dung.
   "She seems to claim if companies start using this method it might help with global warming because some of all the cow dung that causes problems in the atmosphere will start getting used," Abrahams said in an interview.
   
2. Linguistics: Juan Manuel Toro, Josep B. Trobalon and Nuria Sebastian-Galles, of Universitat de Barcelona - for a study showing rats sometimes fail to distinguish between a person speaking Japanese backwards and a person speaking Dutch backwards.
   
3. Peace Prize: The Air Force Wright Laboratory, Dayton, Ohio for instigating research and development on a chemical weapon, the so-called "gay bomb", that "will make enemy soldiers become sexually irresistible to each other".
   
4. Biology: Dr Johanna EMH van Bronswijk of Eindhoven University of Technology, The Netherlands, for their census of all the mites, insects, spiders, pseudoscorpions, crustaceans, bacteria, algae, ferns and fungi that share our beds at night.
   
5. Economics: Kuo Cheng Hsieh, of Taichung, Taiwan, for patenting a device in 2001 that catches bank robbers by dropping a net over them, known as the "net trapping system for capturing a robber immediately".
   The inventor, however, could not be found by Ig Nobel representatives in Taiwan "We had people in Taiwan looking for him. He's vanished. Somebody suggested to us the possibility that maybe the poor man was trapped inside his own machine," Abrahams said.

Enigma 1462

Here is how Mathematica can be used to solve New Scientist Enigma number 1462 in the 29 September 2007 issue.

List of all 3-tuples of side lengths of red, blue, and yellow cubes. The upper limit of each of these is constrained by the known number of digits in the cubes of these side lengths in the red and yellow cases, and by the known number of digits in the square of the side length in the blue case.
(sizelist = Flatten[Table[{r,b,y}, {r,0,9}, {b,0,9}, {y,0,21}], 2]) // Length
2200

Select the cases where red volume is a 3-digit number, and the blue area is a 2-digit number, and the yellow volume is a 4-digit number, and the first digit of the red volume is the same as the first digit of the blue area, and the last digit of the blue area is the same as the last digit of the yellow volume.
sizelist2 = Select[sizelist, 100<=#[[1]]^3<=999 && 10<=#[[2]]^2<=99 && 1000<=#[[3]]^3<=9999 && IntegerDigits[#[[1]]^3][[1]] == IntegerDigits[#[[2]]^2][[1]] && IntegerDigits[#[[2]]^2][[-1]] == IntegerDigits[#[[3]]^3][[-1]]&]
{{5,4,16}, {6,5,15}, {7,6,16}}

Check the digits corresponding to the letters in "nil", "no", and "zero". Only the first case is admissable because it has all digits distinct for all these letters.
Map[
(
{volumered, areablue, volumeyellow} = {#[[1]]^3, #[[2]]^2, #[[3]]^3};
nilnozero = Flatten[Map[IntegerDigits[#]&, {volumered, areablue, volumeyellow}]];
nilozer = Drop[Drop[nilnozero, {-1}], {4}]
)&,
sizelist2
]
{{1,2,5,6,4,0,9}, {2,1,6,5,3,3,7}, {3,4,3,6,4,0,9}}

1. Choose the first case above which thus fixes the values of {n, i, l, o, z, e, r}.
2. Compute the volumes of the 3 coloured cubes.
3. Generate a list of permutations of these volumes, because which pile is which is unknown.
4. List all the possible alternative cases of {t, h, g}, i.e. permutations of one base case.
5. For each of these cases of {t, h, g} compute the total volume "nothing".
6. Generate a list of the numbers of cubes in each of the 3 piles, where the 3rd pile is the one with the unknown number.
{n, i, l, o, z, e, r} = %[[1]]
{volumered, volumeblue, volumeyellow} = {#[[1]]^3, #[[2]]^3, #[[3]]^3}&[sizelist2[[1]]]
volumelist = Permutations[{volumered, volumeblue, volumeyellow}]
thglist = Permutations[Complement[Range[0,9], {n, i, l, o, z, e, r}]]
volumetotallist = Map[FromDigits[{n, o, #[[1]], #[[2]], i, n, #[[3]]}]&, thglist]
numberlist=Append[Map[FromDigits, {{n, o}, {n, o, n, e}}], x]
{1,2,5,6,4,0,9}
{125,64,4096}
{{125,64,4096}, {125,4096,64}, {64,125,4096}, {64,4096,125}, {4096,125,64}, {4096,64,125}}
{{3,7,8}, {3,8,7}, {7,3,8}, {7,8,3}, {8,3,7}, {8,7,3}}
{1637218, 1638217, 1673218, 1678213, 1683217, 1687213}
{16,1610,x}

For each member of the list of volumes (see step 3 above) and each possible total volume (see step 5 above) generate an equation for the number of cubes in the unknown pile, and label each equation with the colour of the unknown pile.
eqns = Flatten[Outer[{Switch[#1[[3]], 125, "r", 64, "b", 4096, "y"], #1.numberlist==#2}&, volumelist, volumetotallist, 1], 1]
{{y,105040+4096 x==1637218}, {y,105040+4096 x==1638217}, {y,105040+4096 x==1673218}, {y,105040+4096 x==1678213}, {y,105040+4096 x==1683217}, {y,105040+4096 x==1687213}, {b,6596560+64 x==1637218}, {b,6596560+64 x==1638217}, {b,6596560+64 x==1673218}, {b,6596560+64 x==1678213}, {b,6596560+64 x==1683217}, {b,6596560+64 x==1687213}, {y,202274+4096 x==1637218}, {y,202274+4096 x==1638217}, {y,202274+4096 x==1673218}, {y,202274+4096 x==1678213}, {y,202274+4096 x==1683217}, {y,202274+4096 x==1687213}, {r,6595584+125 x==1637218}, {r,6595584+125 x==1638217}, {r,6595584+125 x==1673218}, {r,6595584+125 x==1678213}, {r,6595584+125 x==1683217}, {r,6595584+125 x==1687213}, {b,266786+64 x==1637218}, {b,266786+64 x==1638217}, {b,266786+64 x==1673218}, {b,266786+64 x==1678213}, {b,266786+64 x==1683217}, {b,266786+64 x==1687213}, {r,168576+125 x==1637218}, {r,168576+125 x==1638217}, {r,168576+125 x==1673218}, {r,168576+125 x==1678213}, {r,168576+125 x==1683217}, {r,168576+125 x==1687213}}

How many cubes are in the remaining pile and what is their colour? This is the only integer solution that can found to any of the above equations.
Select[Map[{#[[1]], Reduce[#[[2]], x, Integers]}&, eqns], (#[[2]]=!=False&)]
{{b,x==21413}}

Ising model simulation

This posting shows how to use Mathematica (version 6) to create an interactive simulation of a 2-dimensional Ising model, which is based on some work on Boltzmann machines that I did back in the mid 1980s. This nicely illustrates the use of the Manipulate function in Mathematica, which is by far the easiest way of creating interactive graphical demonstrations that I have ever used.

Define a function for randomly initialising the state of the 2-dimensional Ising model.

init[narray_] := array = RandomInteger[{0,1}, {narray,narray}];

Define a function for doing a single Monte Carlo update of the 2-dimensional Ising model:

update[{pe_,pne_,pn_,pnw_}, narray_] :=
Module[
{probpatch={{pnw,pn,pne}, {pe,0.5,pe}, {pne,pn,pnw}}, i, j, di, dj, arraypatch, freq1, freq0},
{i,j} = RandomInteger[{1,narray}, {2}];
arraypatch = Table[array[[Mod[i+di,narray,1],Mod[j+dj,narray,1]]], {di,-1,1}, {dj,-1,1}];
freq1 = Apply[Times, probpatch arraypatch+(1-probpatch)(1-arraypatch), {0,1}];
freq0 = Apply[Times, (1-probpatch)arraypatch+probpatch(1-arraypatch), {0,1}];
array[[i,j]] = If[Random[] < freq1/(freq0+freq1), 1, 0];
];

Define a function for implementing the interactive graphical demonstration.

Manipulate[
If[narray!=Length[array], init[narray], ];
If[u, Do[update[{p2en[[1]],p2nenw[[1]],p2en[[2]],p2nenw[[2]]}, narray], {nupdate}], ];
ArrayPlot[array],
{{p2en, {0.5,0.5}, "probability(\[RightArrow],\[UpArrow])"}, {e, e}, {1-e, 1-e}, ControlPlacement->Left},
{{p2nenw, {0.5,0.5}, "probability(\[UpperRightArrow],\[UpperLeftArrow])"}, {e, e}, {1-e, 1-e}, ControlPlacement->Left},
{{u,"False","update"}, {True,False}, ControlPlacement->Top},
{{narray,32,"size"}, 3, 64, 1, Appearance->"Labeled", ControlPlacement->Top},
FrameLabel->"Ising Model",
Initialization->(e=0.001; array={}; nupdate=50)
]

When Manipulate is evaluated the initial output typically looks like this:

Use the check box at the top of the output window to switch the simulation on and off. The number of Monte Carlo updates that is applied between the display of each "frame" of the simulation is hardwired to be 50.

The initial size of the 2-dimensional array of "spins" is 32 by 32, which is controlled by the 1D slider at the top of the output window. This slider can be moved at any time, even whilst a simulation is being run.

Each Monte Carlo update is implemented by selecting a spin at random, and calculating 8 probability factors due to its interaction with each of its 8 surrounding spins. Each such factor specifies a multiplicative contribution to the relative likelihood that the central spin has the same/different value compared to its neighbouring spin (this is not the most general interaction that is possible). There are 4 independent probability factors corresponding to the following directions from the central spin: east, north, north-east, and north-west. These probability factors are controlled by the two 2D sliders on the left of the output window. The top 2D slider controls the east (left/right slider movement) and north (up/down slider movement) factors, and analogously the bottom 2D slider controls the north-east and north-west factors.

Denote the probability factor for a pair of spins being the same as p, and for being different as 1-p. The central position of each 2D slider corresponds to a probability factor p=1/2 (i.e. the same contribution whether the spins are the same or different), the left and bottom edges of each 2D slider correspond to p=0 (i.e. force the spins to be different), and the right and top edges correspond to p=1 (i.e. force the spins to be the same). These sliders can be moved at any time, even whilst a simulation is being run.

So, using a 1-dimensional example, the spin configuration 0?0 gives a product of probability factors (1-p)^2 for 010 and p^2 for 000, whereas 1?0 gives p(1-p) for 110 and (1-p)p for 100 (i.e. the same in each case). A Monte Carlo update of the central spin selects ? as being 1 or 0 with a relative probability given by the ratio of the corresponding products of probability factors. This ratio is ((1-p)/p)^2 and 1 in the two examples above, respectively. The 2-dimensional case is a straightforward generalisation of the 1-dimensional case.

Each snapshot of the output below shows the typical result of running a randomly initialised simulation for a few seconds with the 2D slider settings as shown in the snapshot.

Note that extreme values are used for the probability factors in each of these simulations, and because the simulations each have a short duration they do not actually reach equilibrium. For instance, the first simulation enforces such strong positive correlations between adjacent spins that it should eventually degenerate to all the spins having the same colour.

Enjoy playing with this Ising model simulation. It is quite interesting to move the 2D sliders to vary the probability factors as the simulation is running, because the speed of the simulation is sufficiently fast that you get an almost real-time response as the Ising model dynamically adjusts its equilibrium state.

High risk, high payoff research

There is an interesting blog posting at Advanced Nanotechnology on The struggle over high risk, high payoff research, which pours scorn on the current low-risk approach to research funding. I have extracted some of the more interesting snippets below.

The blog posting itself makes the following comments:

...the United States science and technology research community has seen a return to a culture which is less likely to pursue high risk/high payoff technology research.

There is a struggle between those who want more High risk, high payoff scientific and technological research and development and those who want only timid, incremental goals who also ridicule even the description of a high payoff technological possibility.

Farber [who] sits on a computer science advisory board at the NSF [says]:

... he has been urging the agency to "take a much more aggressive role in high-risk research." He explains, "Right now, the mechanisms guarantee that low-risk research gets funded. It's always, 'How do you know you can do that when you haven't done it?' A program manager is going to tell you, 'Look, a year from now, I have to write a report that says what this contributed to the country. I can't take a chance that it's not going to contribute to the country.'"

The old head of ARPA Charles Herzfeld says (see here):

...the people that you have to persuade are too busy, don't know enough about the subject and are highly risk-averse ... If the system does not fund thinking about big problems, you think about small problems.

It's all pretty damning stuff. Only slightly tongue-in-cheek, I blame it all on the rise of the use of the spreadsheet as a management and accountancy tool, which makes it far too easy for unimaginative people to become so engrossed with the numbers in their spreadsheets that they overlook the big picture.

Enigma 1461

Here is how Mathematica can be used to solve New Scientist Enigma number 1461 in the 22 September 2007 issue. The solution is optimised for clarity rather than brevity.

Generate a list of all the 24-hour clock squares.
squares = Select[Flatten[Outer[100#1+#2&, Range[0,23], Range[0,59]]], IntegerQ[Sqrt[#]]&]
{0, 1, 4, 9, 16, 25, 36, 49, 100, 121, 144, 225, 256, 324, 400, 441, 529, 625, 729, 841, 900, 1024, 1156, 1225, 1444, 1521, 1600, 1849, 1936, 2025, 2116, 2209, 2304}

Generate all 3-tuples from the list of squares. The convention is that each 3-tuple corresponds to {Tom, Dick, Harry}.
(tuples = Tuples[squares, 3]) // Length
35937

Week 1: filter the list of 3-tuples so that T < D < H and also H = T + D.
(triples1a = Select[tuples, #[[1]]<#[[2]]<#[[3]]&]) // Length
(triples1b = Select[triples1a, #[[3]]==#[[1]]+#[[2]]&]) // Length
5456
5

Use the cyclical relationship between the weeks to obtain the corresponding results for weeks 2 and 3.
triples2b=Map[RotateLeft, triples1b];
triples3b=Map[RotateRight, triples1b];

Extract the case where Harry's square is the same for all weeks, keeping only the value of Harry's square.
harrysquare = Intersection[triples1b, triples2b, triples3b, SameTest->(#1[[3]]==#2[[3]]&)][[1,3]]
400

Extract the 3-tuple of squares {Tom, Dick, Harry} for each week that contain the above value of Harry's square.
tdhsquares = Map[Cases[#, {_,_,harrysquare}][[1]]&, {triples1b, triples2b, triples3b}]
{{144,256,400}, {441,841,400}, {625,225,400}}

Extract Tom's squares.
Transpose[tdhsquares][[1]]
{144, 441, 625}

Drawing on air

I learn from KurzweilAI.net about a report at PhysOrg.com (see here) on a new virtual reality system called "Drawing on Air" that allows artists to create 3-dimensional objects. They say

By putting on a virtual reality mask, holding a stylus in one hand and a tracking device in the other, an artist can draw 3D objects in the air with unprecedented precision. This new system is called “Drawing on Air,” and researchers have designed the interface to be intuitive and provide the necessary control for artists to illustrate complicated artistic, scientific, and medical subjects.

and

... artists can stylize their curves while drawing by dynamically adjusting line thickness and color. Haptic effects enable artists to intuitively adjust line thickness by applying pressure against an imaginary 3D surface, making drawing in the air feel similar to pushing a paintbrush against paper ...

and

... once you have this ability to sketch in the air, there are so many different artistic directions you can go with it ...

I'm quite excited by this development because it gets us closer to having an artistic medium that has the richness of the visualisation "studio" that we each have inside our head. The system is still too expensive for everyday use, so I will use a programmatic approach to 3-dimensional art for the time being (e.g. see the sculpture here), but I already have plans in the pipeline to add a simple "Drawing on Air" type interface to enable the artist to have more intuitive control over the software that creates their artwork, e.g. a standard gamepad can be used to manipulate curved surfaces to form 3-dimensional sculptures, and lots more besides.

Enigma 1460

Here is how Mathematica can be used to solve New Scientist Enigma number 1460 in the 15 September 2007 issue.

This problem is designed to make brute force attack difficult because the set of candidate solutions starts off with (9!)^2=131681894400 elements, i.e. it is the Cartesian product of 2 sets of 9! elements. However, a methodical approach can be used to prune this set down until there is only one case left, i.e. the solution. The main trick is to prune each of the 9! element sets as much as possible before forming the Cartesian product of the residual sets, and keeping track of how many candidate solutions remain at each stage.

Create a list of digits 1,...,9 and a list of all permutations of this.
digits=Range[9]
(digitperms=Permutations[digits])//Length
{1,2,3,4,5,6,7,8,9}
362880

Create a list of the 9 colours and a list of all permutations of this.
colours={"grey","hazel","indigo","jade","khaki","lemon","mauve","navy","orange"};
(colourperms=Permutations[colours])//Length
362880

Without loss of generality, assume that "grey" is at circular position 1, and select all digit permutations where each circular position 2, ... , 9 (i.e. not including the "grey" position) has an odd and even digit on its neighbours.
(digitperms2=Select[digitperms,Apply[And,Map[OddQ[#[[1]]+#[[3]]]&,Most[Partition[#,3,1,{1,1}]]]]&])//Length
2880

Select the colour permutations where "hazel", "indigo", "jade" and "khaki" are in consecutive clockwise circular positions. There is no need to look for wrapped-round cases because "grey" is locked in circular position 1.
(colourperms2=Cases[colourperms,{"grey",___,"hazel","indigo","jade","khaki",___}])//Length
120

Build a Cartesian product list of all remaining digit permutations and colour permutations.
(digitcolourperms=Flatten[Outer[List,digitperms2,colourperms2,1],1])//Length
345600

Select the cases where the "hazel", "indigo" and "jade" digits add up to give the "khaki" digit.
(digitcolourperms2=Select[digitcolourperms,Total[Transpose[Cases[#,{_,"hazel""indigo""jade"}]][[1]]]==Cases[#,{_,"khaki"}][[1,1]]&[Transpose[#]]&])//Length
8640

Select the colour permutations where "lemon", "mauve" and "navy" are in consecutive clockwise circular positions. There is no need to look for wrapped-round cases because "grey" is locked in circular position 1.
(digitcolourperms3=Cases[digitcolourperms2,{_,{"grey",___,"lemon","mauve","navy",___}}])//Length
432

Select the cases where the "lemon" and "mauve" digits add up to give 2 times the "navy" digit.
(digitcolourperms4=Select[digitcolourperms3,Total[Transpose[Cases[#,{_,"lemon""mauve"}]][[1]]]==2Cases[#,{_,"navy"}][[1,1]]&[Transpose[#]]&])//Length
24

Select the cases where the "hazel" digit is the same as the number of times you can find a digit equal to the sum of the digits on its neighbours. This leaves just one possibility which is the required solution.
Select[digitcolourperms4,(x=Transpose[#];Length[Select[Partition[x,3,1,{1,1}],#[[1,1]]+#[[3,1]]==#[[2,1]]&]]==Cases[x,{_,"hazel"}][[1,1]])&][[1]]//Transpose
{{1,grey}, {3,hazel}, {2,indigo}, {4,jade}, {9,khaki}, {5,orange}, {8,lemon}, {6,mauve}, {7,navy}}

Hunting for the fundamental laws of physics

Stephen Wolfram has posted here some interesting news about one of his hobbies - hunting for the fundamental laws of physics, where he outlines how he is both developing and using Mathematica to search for fundamental laws of physics that generate simulated universes which have properties that resemble our own real universe.

The power of his approach lies not only in his use of Mathematica, but more fundamentally in his use of very few axioms to define what the fundamental laws of physics are in the first place. His basic object is a network (i.e. nodes and links-between-nodes), and his basic operation is the mutation of a piece of the network (via the application of a set of rules), which thus allows the network structure to have a dynamical behaviour. In this approach the fundamental laws of physics are determined by the choice of the set of network update rules.

It turns out that various simple consistency criteria cause this approach to give rise to both special relativity and general relativity. That is impressive, starting from a rule-based approach!

One of the challenges is to determine the consequences of a particular choice of network update rules, and to ascertain if they correspond to the behaviour of our known real universe. In general, the network behaviour in response to its update rules can be extremely complicated, and working out what is going on can thus be very difficult and time consuming. The development of Mathematica itself is partly driven by the need to create tools for addressing problems such as this.

Wolfram says that he has not yet found a viable candidate for the fundamental laws of physics using this approach, but that he is hard at work both developing and using Mathematica to achieve this goal. As he says

I certainly think it'll be an interesting - almost metaphysical - moment if we finally have a simple rule which we can tell is our universe. And we'll be able to know that our particular universe is number such-and-such in the enumeration of all possible universes.

I wish him luck in this venture. It would be very impressive if he found that a 3-line Mathematica program was all that was needed to generate the behaviour of our known real universe. Even if he is destined not to discover the fundamental laws of physics using this approach, he will nevertheless have created along the way a very useful toolbox for doing lots of other things, i.e. Mathematica.