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

Improve this question
$\endgroup$
3
  • 1
    $\begingroup$ You might be interested in this question. $\endgroup$
    – J. M.'s persistent exhaustion ♦
    Oct 23, 2012 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, 2012 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, 2012 at 9:17

1 Answer 1

Reset to default
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} *)
Improve this answer
$\endgroup$
6
  • $\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, 2012 at 10:06
  • 1
    $\begingroup$ @arkajad Please see edit, I got rid of the errors. $\endgroup$
    – b.gates.you.know.what
    Oct 23, 2012 at 19:34
  • 1
    $\begingroup$ I find that compilation alone (not to C) is significantly faster. $\endgroup$
    – asim
    Oct 23, 2012 at 21:35
  • $\begingroup$ @asim Good point, thanks. $\endgroup$
    – b.gates.you.know.what
    Oct 24, 2012 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, 2013 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.