Re: Problem with compiled function (is this a bug?)
- To: mathgroup at smc.vnet.net
- Subject: [mg65625] Re: Problem with compiled function (is this a bug?)
- From: Maxim <m.r at inbox.ru>
- Date: Tue, 11 Apr 2006 04:04:41 -0400 (EDT)
- References: <200604080445.AAA25147@smc.vnet.net> <e1ah0k$st$1@smc.vnet.net>
- Sender: owner-wri-mathgroup at wolfram.com
On Sun, 9 Apr 2006 08:35:32 +0000 (UTC), Carl K. Woll <carlw at wolfram.com>
wrote:
> szhorvat at gmail.com wrote:
>> I have a problem with the following compiled function:
>>
>> cpl = Compile[{x,y}, Table[Ceiling[Norm[{i, j} - Ceiling[{
>> x, y}/2]]], {i, 1, x}, {j, 1, y}]]
>>
>> When trying to use this function (for example, as idx = cpl[100,100]),
>> I get the following error:
>>
>> CompiledFunction::"cfte" :
>> Compiled expression 49 Sqrt[2] should be a rank 1 tensor of
>> machine-size integers.
>>
>> CompiledFunction::cfex: External evaluation error at instruction 23;
>> proceeding with uncompiled evaluation.
>>
>> I get the same error even if I try to specify that all variables are
>> Reals:
>>
>> cpl = Compile[{{x, _Real}, {y, _Real}}, Table[Ceiling[Norm[{i, j} -
>> Ceiling[{
>> x, y}/2]]], {i, 1, x}, {j, 1, y}], {{i, _Real}, {j, _Real}}]
>>
>> Is this a bug in Mathematica? If not, what is the right way to compile
>> this function? Does this error occur with earlier versions than 5.2?
>>
>> I'd like to use this function to generate indices for a matrix whose
>> values I want to average radially (ie. I'd like to compute the mean
>> values in the ring around the center of the matrix). Unofrtunately this
>> function takes more than 16 seconds on an 500x500 matrix on my machine
>> which seems unrealistically long for such a simple function. I need to
>> use the function interactively on many datasets, often much larger than
>> 500x500.
>>
>> Szabolcs Horvat
>
> I think the problem is that Norm is a relatively new function (version
> 5) and it hasn't been incorporated into Compile. Witness:
>
> In[1]:=
> cpl=Compile[{x,y},
> Table[Ceiling[Norm[{i,j}-Ceiling[{x,y}/2]]],{i,1,x},{j,1,y}]];
>
> In[2]:=
> cpl[[4]]
>
> Out[2]=
> {{1, 5}, {28, 1, 4}, {28, 0, 5}, {7, -1, 8}, {62, 5, 4, 8, 2, 0},
> {7, 0, 9}, {11, 9, 6}, {4, 18}, {7, 0, 9}, {11, 9, 7}, {4, 14},
> {65, 6, 7, 2, 4}, {65, 0, 1, 3, 1}, {7, 1, 9}, {7, 2, 10},
> {15, 0, 9, 2}, {15, 0, 10, 3}, {92, 260, 3, 0, 2, 3, 0, 3, 3, 0, 4},
> {92, 273, 3, 0, 4, 3, 1, 1, 3, 1, 5}, {91, 49, 3, 1, 5, 2, 1, 1},
> {91, 41, 2, 1, 1, 2, 1, 5}, {51, 4, 5, 4},
> {55, Norm, 2, 1, 4, 2, 1, 5}, {64, 8, 5, 0}, {5, 7, 4, -13},
> {5, 6, 5, -17}, {2}}
>
> You'll see that a Norm is buried amidst the compile code, indicating
> that it wasn't compiled. One workaround is to use (Sqrt[#.#]&) instead
> of Norm, as Dot and Sqrt are compilable:
>
> cpl2 = Compile[{x, y},
> Table[Ceiling[(Sqrt[#.#] &)[{i, j} - Ceiling[{x, y}/2]]], {i, 1, x},
> {j, 1,
> y}]];
>
> In[4]:=
> cpl2[[4]]
>
> Out[4]=
> {{1, 5}, {28, 1, 4}, {28, 0, 5}, {7, 0, 9}, {62, 5, 4, 2, 0},
> {7, 0, 10}, {11, 10, 6}, {4, 21}, {7, 0, 10}, {11, 10, 7}, {4, 17},
> {65, 6, 7, 2, 2}, {65, 0, 1, 3, 1}, {7, 1, 10}, {7, 2, 11},
> {15, 0, 10, 2}, {15, 0, 11, 4}, {92, 260, 3, 0, 2, 3, 0, 4, 3, 0,
> 3}, {92, 273, 3, 0, 3, 3, 1, 1, 3, 1, 3},
> {91, 49, 3, 1, 3, 2, 1, 1}, {91, 41, 2, 1, 1, 2, 1, 3},
> {51, 2, 3, 2}, {89, 2, 2, 10}, {15, 0, 10, 3},
> {91, 55, 3, 0, 3, 3, 0, 2}, {91, 49, 3, 0, 2, 2, 0, 10},
> {63, 9, 10, 0}, {5, 7, 4, -16}, {5, 6, 5, -20}, {2}}
>
> However, it seems to me that you ought to be able to optimize this
> function so that even without Compile it would be much faster. For
> example:
>
> cpl3[x_, y_] :=
> Ceiling[Sqrt[
> N[Outer[Plus, (Range[x] - Ceiling[x/2])^2, (Range[y] -
> Ceiling[y/2])^2]]]]
>
> Some tests:
>
> In[12]:=
> r1=cpl2[1000,1000];//Timing
> r2=cpl3[1000,1000];//Timing
> r1===r2
>
> Out[12]=
> {2.391 Second,Null}
>
> Out[13]=
> {0.125 Second,Null}
>
> Out[14]=
> True
>
> The uncompiled function cpl3 seems to meet your speed requirements.
>
> Carl Woll
> Wolfram Research
>
This is probably version dependent, but in Mathematica 5.2 it's faster to
apply N at the beginning:
In[1]:= cpl3[x_, y_] :=
Ceiling[Sqrt[N[Outer[Plus,
(Range[x] - Ceiling[x/2])^2, (Range[y] - Ceiling[y/2])^2]]]]
In[2]:= cpl3b[x_, y_] :=
Ceiling@ Sqrt@ Outer[Plus,
(N@ Range[x] - Ceiling[x/2])^2, (N@ Range[y] - Ceiling[y/2])^2]
In[3]:= (cpl3[1000, 1000];) // Timing
(cpl3b[1000, 1000];) // Timing
Out[3]= {0.781*Second, Null}
Out[4]= {0.469*Second, Null}
To speed this up further we can exploit the symmetry of the problem:
In[5]:= cpl4[x_, y_] := Drop[#, Boole@ EvenQ@ x, Boole@ EvenQ@ y]&@
Map[Join[Reverse@ Rest@ #, #]&, #, {0, 1}]&@
Ceiling@ Sqrt@ Outer[Plus,
(N@ Range[0, x - Ceiling[x/2]])^2,
(N@ Range[0, y - Ceiling[y/2]])^2]
In[6]:= (cpl4[1000, 1000];) // Timing
Out[6]= {0.141*Second, Null}
For square matrices we can employ something like the Bresenham's
algorithm, where we just keep track of how close we are to the next
circle, and then use the symmetry with respect to the diagonal as well:
In[7]:= cpl5c = Compile[{{n, _Integer}},
Module[{r, d, m, ans},
m = n - Quotient[n + 1, 2] + 1;
ans = Array[0&, {m, m}];
Do[
r = i - 1; d = 0;
Do[
ans[[i, j]] = ans[[j, i]] = r;
If[(d += j + j - 1) > 0,
d -= r + r + 1; r++],
{j, i}],
{i, m}];
ans
]];
cpl5[n_] := If[EvenQ@ n, Drop[#, 1, 1]&, #&]@
Map[Join[Reverse@ Rest@ #, #]&, cpl5c[n], {0, 1}]
In[9]:= (cpl5[1000];) // Timing
Out[9]= {0.078*Second, Null}
Maxim Rytin
m.r at inbox.ru
- References:
- Problem with compiled function (is this a bug?)
- From: szhorvat@gmail.com
- Problem with compiled function (is this a bug?)