# Why doesn't my code compile?

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}]

-
I'm not sure if it is possible to compile such a recursive function. –  Jacob Akkerboom Sep 10 '13 at 18:27
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

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 = DeveloperToPackedArray[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.

-