9
$\begingroup$

My question is related to computing what is called "invariant measure" for a particular well known fractal - the Sierpinski triangle.

We have an array m of four two by two matrices, say

m = {{{1.22951,-0.102459},{0.058548,1.1856}},{{2.34016,2.08943},{-2.31509,1.50439}},
{{1.52091,2.98751},{2.58555,-2.0641}},{{-11857.1,6214.29},{3214.29,-1678.57}}};

and an array a of four vectors, say:

a = {{-0.127561,-0.238217},{-0.322608,0.247723},{-0.601847,-0.165997},{1304.29,-353.571}};

We also have an array c of four coefficients (related to probabilities and determinants of the matrices), for instance:

c = {0.817674, 1.95388, 2.22763, 180.714};

Then we have four functions of a two-dimensional vector x, for i=1,2,3,4 defined as

w[i_,x_] := m[[i]].x + a[[i]]

The recursive formula I am trying to code is

f[ n_, x_ ] := Sum[c[[i]]*f[ n-1, w[i,x] ],{i,1,4}]

with f[ 0 ,x ] = 1 over the unit square 0 < x[[1]] < 1, 0 < x[[2]] < 1, and zero elsewhere:

f[0, x_] := If[x[[1]] > 0 && x[[1]] < 1 && x[[2]] > 0 && x[[2]] < 1, 1, 0]

I would like to plot functions f[ n, x ] for, say, n=1 up to n=10. It takes a lot of recursions, so I would like to find the optimal code that can be compiled for a fast execution. I tried to compile, but all my efforts failed, probably because I do not really understand what can be compiled and how. Without compiling it takes Mathematica 22 seconds to compute just one value::

Timing[f[10, {0.5, 0.5}]]

{22.667, 0.869726}

Can it be done faster with Mathematica? If so, how? Any help will be appreciated

$\endgroup$
  • $\begingroup$ You might be interested in this question. $\endgroup$ – J. M. is away Oct 23 '12 at 7:18
  • $\begingroup$ You can already gain some speed with slight modifications to your definitions : w2[x_] := Dot[m, x] + a ; f2[0, x_] := If[x[[1]] > 0 && x[[1]] < 1 && x[[2]] > 0 && x[[2]] < 1, 1, 0] ; f2[n_, x_] := Dot[c, f2[n - 1, #] & /@ w2[x]]. $\endgroup$ – b.gates.you.know.what Oct 23 '12 at 7:30
  • $\begingroup$ I do not know. Tested again, with a fresh Kernel, got 10 seconds for my original version with w[i,x], 11 seconds with w2[x][[i]]. Thanks for the link, it deals with recursion indeed, but not with compiling. Without compiling, I am afraid, my problem will not fly on Mathematica's wings, I need to cut the time of execution by a factor of at least 100. $\endgroup$ – arkajad Oct 23 '12 at 9:17
9
$\begingroup$

You can very simply compile if you specify explicitly that fc returns a real; this will get rid of the errors. As pointed out by @asim compilation to "C" does not increase speed in this case.

wc = Compile[{{m, _Real, 3}, {a, _Real, 2}, {x, _Real, 1}}, 
  Dot[m, x] + a (*, CompilationTarget -> "C"*)]

fc = Compile[{{m, _Real, 3}, {a, _Real, 2}, {c, _Real, 1}, {n, _Integer}, {x, _Real, 1}},
 If[n == 0, If[x[[1]] > 0 && x[[1]] < 1 && x[[2]] > 0 && x[[2]] < 1, 1, 0], 
      Dot[c, fc[m, a, c, n - 1, #] & /@ wc[m, a, x]]], 
 {{wc[_, _, _], _Real, 2},  {fc[_, _, _, _, _], _Real}}(*, CompilationTarget -> "C"*)]

fc[m, a, c, 10, {0.5, 0.5}] // AbsoluteTiming

(* {2.567999, 0.869726} *)
$\endgroup$
  • $\begingroup$ Thanks. But when I am running: wc = Compile[{{m, _Real, 3}, {a, _Real, 2}, {x, _Real, 1}}, Dot[m, x] + a , CompilationTarget -> "C"], I am getting the message: "A library could not be generated from the compiled function". This is what I was receiving in my original attempts. Some functions would compile to "C", but functions like this one would not. What can be the reason? My CCompilers[] show: {{"Name" -> "Visual Studio", "Compiler" -> CCompilerDriverVisualStudioCompilerVisualStudioCompiler, ... "CompilerName" -> Automatic}} $\endgroup$ – arkajad Oct 23 '12 at 10:06
  • 1
    $\begingroup$ @arkajad Please see edit, I got rid of the errors. $\endgroup$ – b.gates.you.know.what Oct 23 '12 at 19:34
  • 1
    $\begingroup$ I find that compilation alone (not to C) is significantly faster. $\endgroup$ – asim Oct 23 '12 at 21:35
  • $\begingroup$ @asim Good point, thanks. $\endgroup$ – b.gates.you.know.what Oct 24 '12 at 7:14
  • $\begingroup$ This is confusing: Calling CompilePrint@fc (after Needs@"CompiledFunctionTools") shows that wc` and fc still initiate MainEvaluate calls. Is it alright? Is MainEvaluate something necessary that cannot be avoided here? $\endgroup$ – István Zachar Feb 21 '13 at 16:22

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.