{"id":3635,"date":"2011-06-18T16:53:40","date_gmt":"2011-06-18T15:53:40","guid":{"rendered":"http:\/\/www.walkingrandomly.com\/?p=3635"},"modified":"2011-06-18T16:53:40","modified_gmt":"2011-06-18T15:53:40","slug":"making-mathematica-faster-with-compile","status":"publish","type":"post","link":"https:\/\/walkingrandomly.com\/?p=3635","title":{"rendered":"Making Mathematica faster with Compile"},"content":{"rendered":"

Over at Sol Lederman’s fantastic new blog, Playing with Mathematica<\/a>, he shared some code that produced the following figure.<\/p>\n

\"Sol's<\/p>\n

Here’s Sol’s code with an AbsoluteTiming command thrown in.<\/p>\n

f[x_, y_] := Module[{},\r\n  If[\r\n   Sin[Min[x*Sin[y], y*Sin[x]]] >\r\n    Cos[Max[x*Cos[y],\r\n       y*Cos[x]]] + (((2 (x - y)^2 + (x + y - 6)^2)\/40)^3)\/\r\n      6400000 + (12 - x - y)\/30, 1, 0]\r\n  ]\r\n\r\n\u00a0AbsoluteTiming[\r\n\\[Delta] = 0.02;\r\nrange = 11;\r\nxyPoints = Table[{x, y}, {y, 0, range, \\[Delta]}, {x, 0, range, \\[Delta]}];\r\nimage = Map[f @@ # &, xyPoints, {2}];\r\n]\r\nRotate[ArrayPlot[image, ColorRules -> {0 -> White, 1 -> Black}], 135 Degree]<\/pre>\n
This took 8.02 seconds on the laptop I am currently working on (Windows 7 AMD Phenom II N620 Dual core at 2.8Ghz).  Note that I am only measuring how long the calculation itself took and am ignoring the time taken to render the image and define the function<\/em>.<\/p>\n
Compiled functions make Mathematica code go faster<\/strong><\/p>\n
Mathematica<\/a> has a Compile function which does exactly what you’d expect…it produces a compiled version of the function you give it (if it can!).  Sol’s function gave it no problems at all.<\/p>\n
f = Compile[{{x, _Real}, {y, _Real}}, If[\r\n    Sin[Min[x*Sin[y], y*Sin[x]]] >\r\n     Cos[Max[x*Cos[y],\r\n        y*Cos[x]]] + (((2 (x - y)^2 + (x + y - 6)^2)\/40)^3)\/\r\n       6400000 + (12 - x - y)\/30, 1, 0]\r\n   ];\r\n\r\nAbsoluteTiming[\r\n \\[Delta] = 0.02;\r\n range = 11;\r\n xyPoints =\r\n  Table[{x, y}, {y, 0, range, \\[Delta]}, {x, 0, range, \\[Delta]}];\r\n image = Map[f @@ # &, xyPoints, {2}];\r\n ]\r\nRotate[ArrayPlot[image, ColorRules -> {0 -> White, 1 -> Black}],\r\n 135 Degree]<\/pre>\n
This simple change takes computation time down from 8.02 seconds to 1.23 seconds which is a 6.5 times speed up for hardly any extra coding work. Not too shabby!<\/p>\n
Switch to C code to get it even faster<\/strong><\/p>\n
I’m not done yet though!  By default the Compile command produces code for the so-called Mathematica Virtual Machine but recent versions of Mathematica allow us to go even further.<\/p>\n
Install Visual Studio Express 2010 (and the Windows 7.1 SDK if you are running 64bit Windows) and you can ask Mathematica to convert the function to low level C code, compile it and produce a function object linked to the resulting compiled code.  Sounds complicated but is a snap to actually do.  Just add<\/p>\n
CompilationTarget -> \"C\"<\/pre>\n
to the Compile command.<\/p>\n
f = Compile[{{x, _Real}, {y, _Real}},\r\n   If[Sin[Min[x*Sin[y], y*Sin[x]]] >\r\n     Cos[Max[x*Cos[y],\r\n        y*Cos[x]]] + (((2 (x - y)^2 + (x + y - 6)^2)\/40)^3)\/\r\n       6400000 + (12 - x - y)\/30, 1, 0]\r\n   , CompilationTarget -> \"C\"\r\n   ];\r\n\r\nAbsoluteTiming[\\[Delta] = 0.02;\r\n range = 11;\r\n xyPoints =\r\n  Table[{x, y}, {y, 0, range, \\[Delta]}, {x, 0, range, \\[Delta]}];\r\n image = Map[f @@ # &, xyPoints, {2}];]\r\nRotate[ArrayPlot[image, ColorRules -> {0 -> White, 1 -> Black}],\r\n 135 Degree]<\/pre>\n
On my machine this takes calculation time down to 0.89 seconds which is 9 times faster than the original.<\/p>\n
Making the compiled function listable<\/strong><\/p>\n
The current compiled function takes just one x,y pair and returns a result.<\/p>\n
In[8]:= f[1, 2]\r\n\r\nOut[8]= 1<\/pre>\n
It can’t directly accept a list of x values and a list of y values.  For example for the two points (1,2) and (10,20) I’d like to be able to do f[{1, 10}, {2, 20}] and get the results {1,1}.  However what I end up with is an error<\/p>\n
f[{1, 10}, {2, 20}]\r\n\r\nCompiledFunction::cfsa: Argument {1,10} at position 1 should be a machine-size real number. >><\/pre>\n
To fix this I need to make my compiled function listable which is as easy as adding<\/p>\n
RuntimeAttributes -> {Listable}<\/pre>\n
to the function definition.<\/p>\n
f = Compile[{{x, _Real}, {y, _Real}},\r\n   If[Sin[Min[x*Sin[y], y*Sin[x]]] >\r\n     Cos[Max[x*Cos[y],\r\n        y*Cos[x]]] + (((2 (x - y)^2 + (x + y - 6)^2)\/40)^3)\/\r\n       6400000 + (12 - x - y)\/30, 1, 0]\r\n   , CompilationTarget -> \"C\", RuntimeAttributes -> {Listable}\r\n   ];<\/pre>\n
So now I can pass the entire array to this compiled function at once.  No need for Map.<\/p>\n
AbsoluteTiming[\r\n \\[Delta] = 0.02;\r\n range = 11;\r\n xpoints = Table[x, {x, 0, range, \\[Delta]}, {y, 0, range, \\[Delta]}];\r\n ypoints = Table[y, {x, 0, range, \\[Delta]}, {y, 0, range, \\[Delta]}];\r\n image = f[xpoints, ypoints];\r\n ]\r\nRotate[ArrayPlot[image, ColorRules -> {0 -> White, 1 -> Black}],\r\n 135 Degree]<\/pre>\n
On my machine this gets calculation time down to 0.28 seconds, a whopping 28.5 times faster than the original.  Rendering time is becoming much more of an issue than calculation time in fact!<\/p>\n
Parallel anyone?<\/strong><\/p>\n
Simply by adding<\/p>\n
Parallelization -> True<\/pre>\n
to the Compile command I can parallelise the code using threads.  Since I have a dual core machine, this might be a good thing to do.  Let’s take a look<\/p>\n
f = Compile[{{x, _Real}, {y, _Real}},\r\n   If[\r\n    Sin[Min[x*Sin[y], y*Sin[x]]] >\r\n     Cos[Max[x*Cos[y],\r\n        y*Cos[x]]] + (((2 (x - y)^2 + (x + y - 6)^2)\/40)^3)\/\r\n       6400000 + (12 - x - y)\/30, 1, 0]\r\n   , RuntimeAttributes -> {Listable}, CompilationTarget -> \"C\",\r\n   Parallelization -> True];\r\n\r\nAbsoluteTiming[\r\n \\[Delta] = 0.02;\r\n range = 11;\r\n xpoints = Table[x, {x, 0, range, \\[Delta]}, {y, 0, range, \\[Delta]}];\r\n ypoints = Table[y, {x, 0, range, \\[Delta]}, {y, 0, range, \\[Delta]}];\r\n image = f[xpoints, ypoints];\r\n ]\r\nRotate[ArrayPlot[image, ColorRules -> {0 -> White, 1 -> Black}],\r\n 135 Degree]<\/pre>\n
The first time I ran this it was SLOWER than the non-threaded version coming in at 0.33 seconds.  Subsequent runs varied and occasionally got as low as 0.244 seconds which is only a few hundredths of a second faster than the original.<\/p>\n
If I make the problem bigger, however, by decreasing the size of Delta then we start to see the benefit of parallelisation.<\/p>\n
AbsoluteTiming[\r\n \\[Delta] = 0.01;\r\n range = 11;\r\n xpoints = Table[x, {x, 0, range, \\[Delta]}, {y, 0, range, \\[Delta]}];\r\n ypoints = Table[y, {x, 0, range, \\[Delta]}, {y, 0, range, \\[Delta]}];\r\n image = f[xpoints, ypoints];\r\n ]\r\nRotate[ArrayPlot[image, ColorRules -> {0 -> White, 1 -> Black}],\r\n 135 Degree]<\/pre>\n
The above calculation (sans rendering) took 0.988 seconds using a parallelised version of f and 1.24 seconds using a serial version.  Rendering took significantly longer!  As a comparison lets put a Delta of 0.01 in the original code:<\/p>\n
f[x_, y_] := Module[{},\r\n  If[\r\n   Sin[Min[x*Sin[y], y*Sin[x]]] >\r\n    Cos[Max[x*Cos[y],\r\n       y*Cos[x]]] + (((2 (x - y)^2 + (x + y - 6)^2)\/40)^3)\/\r\n      6400000 + (12 - x - y)\/30, 1, 0]\r\n  ]\r\n\r\n AbsoluteTiming[\r\n\\[Delta] = 0.01;\r\nrange = 11;\r\nxyPoints = Table[{x, y}, {y, 0, range, \\[Delta]}, {x, 0, range, \\[Delta]}];\r\nimage = Map[f @@ # &, xyPoints, {2}];\r\n]\r\nRotate[ArrayPlot[image, ColorRules -> {0 -> White, 1 -> Black}], 135 Degree]<\/pre>\n
The calculation time (again, ignoring rendering time) took 32.56 seconds and so our C-compiled, parallel version is almost 33 times faster!<\/p>\n
Summary<\/strong><\/p>\n