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.
Clear["`*"];
cf = Compile[{{A, _Real, 2}, {n, _Integer}},
   With[{a = A[[1]], b = A[[2]], c = A[[3]]},
    If[n < 1, {{a, b, c}},
      Join[cf[{a, (a + b)/2, (a + c)/2}, n - 1],
       cf[{(b + a)/2, b, (b + c)/2}, n - 1],
       cf[{(c + a)/2, (c + b)/2, c}, n - 1]
       ]
     ];
    ]
   ];

cf[{{0, 0}, {1, 0}, {.5, .8}} // N, 3]

This code I wrote can't be compiled, Mathematica returns

CompiledFunction::cfex: Could not complete external evaluation at instruction 19; proceeding with uncompiled evaluation. >>

How can I correctly compile it?

The uncompiled version:

Clear["`*"];
f[{a_, b_, c_}, n_] :=
  If[n < 1, {{a, b, c}},
   Join[f[{a, (a + c)/2, (a + b)/2}, n - 1],
    f[{(b + a)/2, b, (b + c)/2}, n - 1],
    f[{(c + a)/2, (c + b)/2, c}, n - 1]]
   ];

n = 3;
data = f[{{0, 0}, {1, 0}, {.5, .8}} // N, n];
ListLinePlot[data /. {a_, b_, c_} :> {a, b, c, a}]
share|improve this question
1  
I'm not sure if it is possible to compile such a recursive function. –  Jacob Akkerboom Sep 10 '13 at 18:27
4  
There are several problems with your function: think about what is the type of the return value. What rank tensor is it. I can't really figure out what you're trying to achieve so I can't correct it for you, but: 1. as it is, your function doesn't return anything because of the ; 2. if we fix that, it returns a rank 3 tensor which is fed back as input after a transformation. The input is expected to be a rank 2 tensor however. I'd suggest to first write an uncompiled version and fix the bugs. –  Szabolcs Sep 10 '13 at 18:31
    
@JacobAkkerboom I thought for a moment I could game that restriction, but it seems I was wrong. So I agree, that's the main restriction here. What I pointed out is secondary. –  Szabolcs Sep 10 '13 at 18:37
    
Here it's shown how to create a recursive compiled function, but it relies in MainEvaluate so I'm not sure how much compilation would help with efficiency, unless the non-recursive part is substantial. –  Szabolcs Sep 10 '13 at 18:39
    
Related: (13504), (31171) –  Mr.Wizard Sep 10 '13 at 21:13

1 Answer 1

up vote 3 down vote accepted

Here is a compiled version. First note, however, that your compiled version and the top level version are not the same.

ClearAll[cf]
cf = Compile[{{A, _Real, 2}, {n, _Integer}}, 
   Block[{a = A[[1]], b = A[[2]], c = A[[3]]},
    If[n < 1, A,
     Join[cf[{a, (a + c)/2, (a + b)/2}, n - 1], 
      cf[{(b + a)/2, b, (b + c)/2}, n - 1], 
      cf[{(c + a)/2, (c + b)/2, c}, n - 1]]
     ]
    ]
   ];

This is called via (see further down):

dc = Partition[cf[input, n], 3];

A slightly different version that does not need an additional Partition

ClearAll[cf2]
cf2 = Compile[{{A, _Real, 2}, {n, _Integer}}, 
   Block[{a = A[[1]], b = A[[2]], c = A[[3]]},
    If[n < 1, {A},
     Join[cf2[{a, (a + c)/2, (a + b)/2}, n - 1], 
      cf2[{(b + a)/2, b, (b + c)/2}, n - 1], 
      cf2[{(c + a)/2, (c + b)/2, c}, n - 1]]
     ]
    ]
   , {{cf2[_, _], _Real, 3}}
   ];

For a comparison:

n = 10;
input = Developer`ToPackedArray[N[{{0, 0}, {1, 0}, {.5, .8}}]];
AbsoluteTiming[dc = Partition[cf[input, n], 3];]
AbsoluteTiming[dc2 = cf2[input, n];]
AbsoluteTiming[data = f[input, n];]
data === dc === dc2

(*
0.18
0.18
0.71
True
*)

Also, a better way to rearrange the data is:

data[[All, {1, 2, 3, 1}]]

which does not unpack.

ListLinePlot[data[[All, {1, 2, 3, 1}]]]

If you look at the compiled code with CompilePrint you will note that there are calls to CopyTensor - with some clever computed array (by direct assignment of the recursive cf calls into an array in stead of using Join) it may be possible to avoid those and save computational expense.

Hope this helps.

share|improve this answer

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.