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.

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]
share|improve this question
2  
Just put the //Abs//Log into the Compile –  Simon Woods Nov 18 '13 at 13:56

3 Answers 3

up vote 11 down vote accepted

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.

share|improve this answer
    
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 –  Stefan Nov 18 '13 at 16:35

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.

share|improve this answer

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.

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.