9
$\begingroup$

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?

enter image description here

 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]
$\endgroup$
2

3 Answers 3

11
$\begingroup$

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

enter image description here

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

$\endgroup$
1
  • $\begingroup$ I think this is a pretty good answer, since it shows that you don't have to pull always the big gun (Compile), if you think about proper data structures. Well done, indeed. +1 $\endgroup$
    – Stefan
    Nov 18, 2013 at 16:35
8
$\begingroup$

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.

$\endgroup$
1
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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