Compiling ColorFunction for faster complex phase-amplitude plots

Question

I have gotten into the habit of plotting complex-valued functions on a plane using a color-function along the lines of

Function[{z}, Hue[Arg[z]/(2 π), 1, Abs[z]]].

For example, the following code generates an interesting image of the sum of 5 spherical wave sources (without 1/r decay included):

i = 130;
ω = 2.0 π i/50;
n = 5;
L1 = 3.0 - δ/2;
L2 = 4.0 - δ/2;
δ = 0.01;
f = Compile[{{x, _Real}, {y, _Real}}, 
Evaluate[Sum[Exp[ω I Sqrt[(x - 1.0 Cos[θ])^2 + (y - 
      1.0 Sin[θ])^2]], {θ, 2 π/n, 2 π, 
  2 π/n}]]];
ArrayPlot[
 Table[f[x, y], {x, -L1, L1, δ}, {y, -L2, L2, δ}], 
 Frame -> False, 
 ColorFunction -> Function[{z}, Hue[Arg[z]/(2 π), 1, Abs[z/2]]], 
 ColorFunctionScaling -> False, PixelConstrained -> {1, 1}]

I have two questions. The second is a little more important than the first.

1. Exchanging ArrayPlot with ListDensityPlot in the above example either simply yields a blank image or flat-out crashes my MathKernel. I believe this is because ListDensityPlot sees the complex-valued array, freaks out because they're not real, and spits out a blank plot even though it has a user-supplied color-function for handling the complex entries of the array. Is there any way of getting ListDensityPlot to properly use this color-function?

2. The example I gave above is pretty slow. The construction of the array alone takes 1.26 seconds on my machine, whereas the displaying of the code cited above takes 10.27 seconds (it's just an 800 by 600 pixel plot). Replacing the phase-amplitude ColorFunction I used above with the built-in function GrayLevel (and applying Abs to the array) gives a timing of 2.43 seconds, and using the default ArrayPlot scheme by removing the references to ColorFunction altogether executes in 1.52 seconds. Subtracting the time to build the array itself, the phase-amplitude color function generates a plot in 9.01 seconds, GrayLevel executes in 1.17 seconds, and the default scheme executes in 0.26 seconds. I am attempting to generate a movie, and 10 seconds to generate each frame is a bit on the slow side. Is there a way to compile this color function to make it run as fast as default, or perhaps as fast as GrayLevel?

For reference purposes, I am running on a 2008 Penryn Dual-Core MacBook and Mathematica 9.0.

Answer 1

Another problem is that the Table[..., {x, -L1, L1, δ}, {y, -L2, L2, δ}] produces unpacked array.

f = With[{fun = 
     Evaluate[
       Sum[Exp[ω I Sqrt[(#1 - 1.0 Cos[θ])^2 + (#2 - 
               1.0 Sin[θ])^2]], {θ, 2 Pi/n, 2 Pi, 
         2 Pi/n}]] &}, Compile[{{x, _Real}, {y, _Real}},
    {Mod[Arg[#]/(2.0 Pi), 1], 1.0, Abs[#/2]} &@fun[x, y], 
    RuntimeAttributes -> {Listable}, RuntimeOptions -> "Speed", CompilationTarget -> "C"]];

Image[f[#[[All, All, 1]], #[[All, All, 2]]] &@
  Outer[List, Range[-L1, L1, δ], Range[-L2, L2, δ]], 
   ColorSpace -> "HSB"]

This solution is 4-5 times faster than halirutan's approach for me.

Answer 2

Regarding your 1. question

A density plot is clearly not recommended for your problem. Firstly, a density is AFAIK per definition not complex, but let's ignore this for a moment. The real neck-breaker here is, that ListDensisityPlot interpolates the values if you don't turn it explicitly off. And even if you turn it off, a ListDensityPlot will create a Graphics with polygons. This will make the whole thing very slow. My recommendation is: if you have a fixed equidistant grid (as opposed to: you need a variably refined grid) you should always try to use Image, because it has incredible performance advantages.

Regarding your 2. question:

I would do the following: Compile your complex sum and the coloring into one big function, which is parallelized and can be called directly on your large matrix. For this I made the argument to your function complex and transferred the hue value calculation directly to the end:

i = 130;
ω = 2.0 π i/50;
n = 5;
L1 = 3.0 - δ/2;
L2 = 4.0 - δ/2;
δ = 0.01;
f = Compile[{{z, _Complex, 0}}, Evaluate[
    With[{res =
       Sum[
        Exp[ω I Sqrt[(Re[z] - 1.0 Cos[θ])^2 + (Im[z] - 
               1.0 Sin[θ])^2]], {θ, 2 π/n, 2 π, 
         2 π/n}]},
     {Arg[res]/(2 π) + 1/2, 1, Abs[res/2]}
     ]
    ],
CompilationTarget -> "C", Parallelization -> True, RuntimeAttributes -> {Listable}];

And now you make just

Image[f@Table[x + I y, {x, -L1, L1, δ}, {y, -L2, L2, δ}], ColorSpace -> "HSB"]

Which runs in 0.383 seconds here for everything except the compilation step. The advantage over ssch's approach (which needs 0.63 seconds on my machine) is, that f@Table[..] uses low-level parallelization, which is very fast if applied correctly.

The fact that Table does not produce packed arrays (which are used in compiled C functions) is already known for some time. I remember that we already have a Q&A about this. Therefore, one way to improve this further is to use e.g. Outer as already pointed out by ybeltukov. Although, the call with Outer and my complex function is (IMO) easier to read than the array access of ybeltukov

Image[f@Outer[Complex, Range[-L1, L1, δ], Range[-L2, L2, δ]], ColorSpace -> "HSB"]
(* compared to 
Image[f[#[[All, All, 1]], #[[All, All, 2]]] &@
  Outer[List, Range[-L1, L1, δ], Range[-L2, L2, δ]], 
   ColorSpace -> "HSB"]
*)

it comes with a drawback: Using a compiled function which works on complex numbers directly (instead of x and y as real and imaginary part) seems always to be slower (even if I adapt my function to look very close to ybeltukov's implementation).

Answer 3

Here is an approach to create an Image which gets its speed from using the listability of Arg and Abs.

EDIT Updated with Outer from @ybeltukov and Parallelization from @halirutan. It is now has the same speed as ybeltukov :)

This method has the advantage that the color function doesn't have to be compiled into main function:

(* Same as in originial question with new Compile options *)
f = Compile[{{x, _Real}, {y, _Real}}, 
   Evaluate[
     Sum[Exp[ω I Sqrt[(x - 1.0 Cos[θ])^2 + (y - 1.0 Sin[θ])^2]], {θ, 2 π/n, 2 π, 2 π/n}]],
   RuntimeAttributes -> {Listable},
   CompilationTarget -> "C",
   Parallelization -> True];

(* Slower, used in halirutans comparision:
   grid = Table[f[x, y], {x, -L1, L1, δ}, {y, -L2, L2, δ}];*)

(* Improved with Outer from ybeltukovs answer *)
grid = f[#[[All, All, 1]],
         #[[All, All, 2]]] &@Outer[List, Range[-L1, L1, δ], Range[-L2, L2, δ]];
imgdata={
  Arg[grid]/(2.Pi)+0.5,
  ConstantArray[1, Dimensions[grid]],
  Abs[grid]/2};

Image[imgdata, ColorSpace -> "HSB", Interleaving -> False]
(* 
 Table: 1.0s not used anymore
 Outer: 0.18s
 imgdata: 0.05
 Image: 0.002s *)

Answer 4

Clear["`*"];
cf = Compile[{},
   Module[{i, ω, n, L1, L2, δ, dat},
    i = 130;
    ω = 2.0 π i/50;
    δ = 0.01;
    n = 5;
    L1 = 3.0 - δ/2;
    L2 = 4.0 - δ/2;
    dat = With[{y = Range[-L2, L2, δ]},
      Table[
       Sum[Exp[ω I Sqrt[(x - 1.0 Cos[θ])^2 + (y - 
               1.0 Sin[θ])^2]], {θ, 2 Pi/n, 2 Pi, 
         2 Pi/n}],
       {x, -L1, L1, δ}]];

    {Mod[Arg[#]/(2.0 Pi), 1], ConstantArray[1.0, Dimensions@#], 
        Abs[#/2]} &@dat~Transpose~{3, 1, 2}

    ]
   ];
data = cf[]; // AbsoluteTiming
Image[data, ColorSpace -> "HSB"]

Answer 5

Here's one with okay performance without using Compile. I create lists of x and y, then construct the sum, keeping x and y as 1D lists until it's necessary to combine them with Outer. This saves some calculation.

Module[{x, y, q, z},
  x = Range[-L1, L1, δ];
  y = Range[-L2, L2, δ];
  q = Sum[Exp[ω I Sqrt[
       Outer[Plus, (x - 1.0 Cos[θ])^2, (y - 1.0 Sin[θ])^2]
       ]], {θ, 2 π/n, 2 π, 2 π/n}];
  z = ConstantArray[1, {Length[x], Length[y]}];
  Image[{Arg[q]/(2 π) + 0.5, z, Abs[q/2]}, ColorSpace -> "HSB", Interleaving -> False]]