Compiling the VoigtDistribution PDF

Question

According to List of compilable functions, Erf and Erfc are compilable functions.

However, I want to make a compiled version of the PDF of a VoigtDistribution to use in a NonlinearModelFit, and it doesn't seem that the Erfc of a complex value will compile:

funcReal = 
 Compile[{{x, _Real}}, Erfc[x I], CompilationTarget -> "C", 
  RuntimeOptions -> "Speed"]
funcComplex = 
 Compile[{{x, _Complex}}, Erfc[x I], CompilationTarget -> "C", 
  RuntimeOptions -> "Speed"]

Needs["CompiledFunctionTools`"];
CompilePrint[funcReal] (*same as funcComplex*)

    1 argument
    2 Real registers
    4 Complex registers
    Underflow checking off
    Overflow checking off
    Integer overflow checking off
    RuntimeAttributes -> {}

    R0 = A1
    C0 = 0. + 1. I
    R1 = 0.
    Result = C3

1   C1 = R0 + R1 I
2   C1 = C1 * C0
3   C2 = R0 + R1 I
4   C2 = C2 * C0
5   C3 = MainEvaluate[ Hold[Erfc][ C2]]
6   Return

Note the call to MainEvaluate:

Erfc[I] // N
funcReal[1]
funcComplex[1]

1. - 1.65043 I

1. - 1.65043 I

1. - 1.65043 I

All the functions work, but because of the MainEvaluate, they offer no performance benefit. How can I compile this function? Is this possible? Is there an alternative formula I could use?

Removing the CompilationTarget doesn't solve the problem either.

They compile for real values, but do not appear to compile for complex arguments. — asim, Feb 20 '13 at 17:56
Look at e.g. Plot3D[Arg@Erf[x + I y], {x, -5, 5}, {y, -5, 5}]. This is not a function you would probably wish to try to find an alternative formula for, even approximately. — Oleksandr R., Feb 20 '13 at 18:25
On the other hand, regarding the PDF of the Voigt distribution: http://dx.doi.org/10.1016/0368-2048(94)02189-7 — Oleksandr R., Feb 20 '13 at 18:46

score 11 · Answer 1 · answered Nov 13 '15 at 14:51

Here is a compiled implementation of the Voigt profile function, based on an approximation derived by Chiarella and Reichel and improved by Abrarov, Quine and Jagpal:

voigt = With[{n = 24, τ = 12},
             With[{d = N[Range[n] π/τ], b = N[Exp[-(Range[n] π/τ)^2]],
                   s = N[PadRight[{}, n, {-1, 1}]],
                   sq = N[Sqrt[2]], sp = N[Sqrt[2 π]]},
                  Compile[{{δ, _Real}, {σ, _Real}, {x, _Real}},
                          Module[{z = (x + I δ)/(σ sq), e}, e = Exp[I τ z];
                                 Re[(I (1 - e)/(τ z) + (2 I z/τ)
                                    b.((e s - 1)/((d + z) (d - z))))]/(σ sp)],
                          RuntimeAttributes -> {Listable}]]];

The compiled function performs remarkably well on my box:

vp = PDF[VoigtDistribution[1, 3/2]];
AbsoluteTiming[v1 = vp[N[Range[-7, 7, 1/50]]];]
   {3.70813, Null}

AbsoluteTiming[v2 = voigt[1, 3/2, Range[-7, 7, 1/50]];]
   {0.053813, Null}

Norm[v1 - v2, ∞]
   1.38778*10^-16

Further speed could probably be squeezed out by explicitly expanding out the complex arithmetic into its real and imaginary components, but the resulting code will look comparatively unwieldy.

You could consider submitting this to the new function repository: https://resources.wolframcloud.com/FunctionRepository/ It would be a great addition! — Sjoerd Smit, Mar 15 '19 at 16:44

score 8 · Accepted Answer · edited Nov 17 '15 at 18:01

From the documentation, we have a pseudo-Voigt Distribution that might be used as an approximation. This might be useful as a basis for making a compilable function.

PseudoVoigtDistribution[δ_, σ_] := 
 Block[{g = (δ^5 + σ^5 + 2.69296 σ^4 δ + 2.42843 σ^3 δ^2 + 
        4.47163 σ^2 δ^3 + 0.07842 σ δ^4)^(1/5), η},
       η = δ/g;
       η = η*(1.36603 - 0.47719 η + 0.11116 η^2);
       MixtureDistribution[{1 - η, η}, {NormalDistribution[0, g], 
                           CauchyDistribution[0, g]}]
       ]

Compiling the VoigtDistribution PDF

2 Answers2

Linked