Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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

share|improve this question
    
You might be interested in this question. –  J. M. Oct 23 '12 at 7:18
    
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]]. –  b.gatessucks Oct 23 '12 at 7:30
    
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. –  arkajad Oct 23 '12 at 9:17
add comment

1 Answer 1

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} *)
share|improve this answer
    
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" -> CCompilerDriver`VisualStudioCompiler`VisualStudioCompiler, ... "CompilerName" -> Automatic}} –  arkajad Oct 23 '12 at 10:06
1  
@arkajad Please see edit, I got rid of the errors. –  b.gatessucks Oct 23 '12 at 19:34
1  
I find that compilation alone (not to C) is significantly faster. –  asim Oct 23 '12 at 21:35
    
@asim Good point, thanks. –  b.gatessucks Oct 24 '12 at 7:14
    
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? –  István Zachar Feb 21 '13 at 16:22
show 1 more comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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