What is wrong with my expression for compiling?

Question

Following code was translated from a snippet matlab code. Without compiling it worked slowly. I want to compile it, but am getting a numerical error. How do I compile this in a correct way?

 data = Compile[{}, 
      Table[(1/49 1. x^2 Sqrt[Abs[Abs[x] - 3]/(Abs[x] - 3)] + 
          1/9 y^2 Sqrt[ Abs[y + (3 Sqrt[33])/7]/(y + (3 Sqrt[33])/7)] - 
          1) (Abs[x]/2 - 1/112 (3 Sqrt[33] - 7) x^2 - 3 + 
          Sqrt[1 - (Abs[Abs[x] - 2] - 1)^2] - 
          y) (9 Sqrt[Abs[(Abs[x] - 1) (Abs[x] - 0.75)]/((1 - Abs[x]) (Abs[x] - 
            0.75))] - 8 Abs[x] - y) (3 Abs[x] + 
          0.75 Sqrt[ Abs[(Abs[x] - 0.75) (Abs[x] - 0.5)]/((0.75 -  Abs[x]) (Abs[x] -
           0.5))] -  y) (2.25 Sqrt[ Abs[(x - 0.5) (x + 0.5)]/((0.5 - x)
          (x + 0.5))] -  y) ((6 Sqrt[10])/ 7 + (1.5 - 0.5 Abs[x]) Sqrt[
            Abs[Abs[x] - 1]/(Abs[x] - 1)] - 
         3/7 Sqrt[10] Sqrt[4 - (Abs[x] - 1)^2] - y),
         {y, 4.5, -4.5, -0.041}, {x, -7, 7, 0.041}]
      ][] // Abs // Log;
ArrayPlot[data]

Answer 1

Often compilation isn't necessary if you use packed arrays. Here is my attempt with some simplifications:

batman[dx_: 0.041, dy_: 0.041, x_: 7, y_: 4.5] := 
  ArrayPad[#, {{0}, {Dimensions[#][[2]] - 1, 0}}, 
         Padding -> "Reflected"] &@
       Evaluate@
        Log@Abs@N[-(1/
             2765952) (8 #1 - 
              18 Sqrt[Abs[(-1 + #1) (-(3/4) + #1)]] Sqrt[(
               1 + 0.0 I)/(-3 + 7 #1 - 4 #1^2)] + #2) (-9 Sqrt[
               Abs[1/2 - #1]] Sqrt[(1 + 0.0 I)/(2 - 4 #1)] + 
              2 #2) (-6 #1 - 
              3 Sqrt[2] Sqrt[Abs[3/8 - (5 #1)/4 + #1^2]] Sqrt[(
               1 + 0.0 I)/(-3 + 10 #1 - 8 #1^2)] + 
              2 #2) (-112 Sqrt[2 + 0.0 I - Abs[-2 + #1]] Sqrt[
               Abs[-2 + #1]] - 56 #1 + (-7 + 3 Sqrt[33]) #1^2 + 
              112 (3 + #2)) (7 Sqrt[Abs[-1 + #1]] (-3 + #1) Sqrt[(
               1 + 0.0 I)/(-1 + #1)] + 
              2 (3 Sqrt[10] (-2 + Sqrt[3 + 0.0 I + 2 #1 - #1^2]) + 
                 7 #2)) (9 Sqrt[Abs[-3 + #1]] Sqrt[(
               1 + 0.0 I)/(-3 + #1)] #1^2 + 
              49 (-9 + 
                 Sqrt[7] Sqrt[Abs[(3 Sqrt[33])/7 + #2]] #2^2 Sqrt[(
                  1 + 0.0 I)/(3 Sqrt[33] + 7 #2)]))] & @@ 
     Transpose[#, {3, 2, 1}] &@
   Outer[List, Range[0, x, dx], Range[y, -y, -dy]];

Its timing is comparable to Peltio's compiled solution (even a bit faster because I calculate only one of the symmetric parts)

data = batman[]; // AbsoluteTiming

{0.078602, Null}

It produces packed array

data // Developer`PackedArrayQ

True

Now one can obtain hi-res image

Note: I add 0.0 I in square roots because Sqrt of a real array with negative signs can unpack it.

Answer 2

Well, I guess you should give a couple of arguments to Compile, but first you could rewrite your expression in a slightly simplified version:

batsign[x_, y_] = FullSimplify[(1/
    49 1. x^2 Sqrt[Abs[Abs[x] - 
      3]/(Abs[x] - 3)] + 1/9 y^2 Sqrt[Abs[
                y + (3 Sqrt[33])/7]/(y + (3 Sqrt[33])/7)] - 1) (Abs[x]/2 - 1/
            112 (3 Sqrt[33] - 7) x^2 - 3 + Sqrt[
                    1 - (Abs[Abs[x] - 2] - 1)^2] - y) (9 Sqrt[Abs[(Abs[
                      x] - 1) (Abs[x] - 0.75)]/((1 - Abs[
            x]) (Abs[x] - 0.75))] - 8 Abs[
            x] - y) (3 Abs[x] + 0.75 Sqrt[Abs[(Abs[x] - 0.75) (Abs[x] - 
            0.5)]/((0.75 - Abs[
              x]) (Abs[x] - 0.5))] - y) (2.25 Sqrt[Abs[(x - 0.5) (
                  x + 0.5)]/((0.5 - x) (x + 0.5))] - y) ((
                    6 Sqrt[10])/7 + (1.5 - 0.5 Abs[x]) Sqrt[Abs[
                          Abs[x] - 1]/(Abs[x] - 1)] - 
                          3/7 Sqrt[10] Sqrt[4 - (Abs[x] - 1)^2] - y)];

This might reveal some common factors to speed up the compiled version a hair more. Just use the returned expression (remove the semicolon at the end to show it) instead of batsign[x,y] below, if you do not want to simplify it every time (let me clarify: this means that you won't have to account for the time taken by FullSimplify to do its job). This will compile

batsignC[x_, y_] = Compile[{x, y}, Log[Abs[batsign[x, y]]]];

data = Table[batsignC[x, y],
        {y, 4.5, -4.5, -0.041}, {x, -7, 7, 0.041}
    ];

Then you can plot it

ArrayPlot[data]

Is it faster? Use Timing to get the time required to compute data in both cases and compute the ratio. This should give some sort of machine independent measure of how much faster the compiled code is.

Answer 3

Thanks for @Simon Woods, this one would be faster than that Matlab code:

cf = Compile[{}, 
   With[{y = Range[4.5, -4.5, -0.02 + 10.^-6]}, 
    Table[((x/7)^2 Sqrt[Abs[Abs[x] - 3]/(Abs[x] - 3)] + (y/3)^2 Sqrt[
            Abs[y + (3 Sqrt[33])/7]/(y + (3 Sqrt[33])/7)] - 
          1) (Abs[x]/2 - 1/112 (3 Sqrt[33] - 7) x^2 - 3 + 
          Sqrt[1 - (Abs[Abs[x] - 2] - 1)^2] - y) (9 Sqrt[
            Abs[(Abs[x] - 1) (Abs[x] - 0.75)]/((1 - Abs[x]) (Abs[x] - 
                 0.75))] - 8 Abs[x] - y) (3 Abs[x] + 
          0.75 Sqrt[ Abs[(Abs[x] - 0.75) (Abs[x] - 0.5)]/((0.75 - 
                 Abs[x]) (Abs[x] - 0.5))] -  y) (2.25 Sqrt[
            Abs[(x - 0.5) (x + 0.5)]/((0.5 - x) (x + 0.5))] - 
          y) ((6 Sqrt[10])/7 + (1.5 - 0.5 Abs[x]) Sqrt[Abs[Abs[x] - 1]/
       (Abs[x] - 1)] -  3/7 Sqrt[10] Sqrt[4 - (Abs[x] - 1)^2] - y),
 {x, -7, 7, 0.02 + 10.^-6}] // Log // Abs]];
data = cf[] // Transpose; // AbsoluteTiming
Colorize[Image@Rescale@data, ColorFunction -> "DeepSeaColors"]

{0.192011, Null}

Matlab version:

tic
[x, y] = meshgrid(-7: 0.02+1e-6 :7, 4.5: -0.02+1e-6 : -4.5);
batman = (((x/7).^2.*sqrt(abs(abs(x)-3)./(abs(x)-3))+((y/3).^2) .*  ...
           sqrt(abs(y+(3*sqrt(33)/7))./(y+(3*sqrt(33)/7))))-1) .* ...
         (abs(x/2)-((3*sqrt(33)-7)/112).*(x.^2)-3+sqrt(1-(abs(abs(x)-2)-1).^2) - y) .* ...
         (9*sqrt(abs((abs(x)-1).*(abs(x)-0.75))./((1-abs(x)) .* (abs(x)-0.75)))-8*abs(x)-y) .* ...
         (3*abs(x)+0.75*sqrt(abs((abs(x)-0.75).*(abs(x)-0.5))./((0.75-abs(x)).*(abs(x)-0.5)))-y) .* ...
         (2.25*sqrt(abs((x-0.5).*(x+0.5))./((0.5-x).*(0.5+x)))-y) .* ...
         (((6*sqrt(10))/7)+(1.5-0.5*abs(x)) .* sqrt(abs(abs(x)-1)./(abs(x)-1))-((6*sqrt(10))/14).*sqrt(4-(abs(x)-1).^2)-y);
toc
imshow( log(abs(batman)), [])

Elapsed time is 0.393433 seconds.