Compiling the VoigtDistribution PDF

13

1

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.

s0rce

Posted 2013-02-20T17:17:45.750

Reputation: 9 292

3They compile for real values, but do not appear to compile for complex arguments. – asim – 2013-02-20T17:56:50.377

1Look 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. – 2013-02-20T18:25:59.387

1

On the other hand, regarding the PDF of the Voigt distribution: http://dx.doi.org/10.1016/0368-2048(94)02189-7

– Oleksandr R. – 2013-02-20T18:46:33.120

Answers

7

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]}]
       ]

asim

Posted 2013-02-20T17:17:45.750

Reputation: 1 855

9

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.

J. M.'s ennui

Posted 2013-02-20T17:17:45.750

Reputation: 115 520

1

You could consider submitting this to the new function repository: https://resources.wolframcloud.com/FunctionRepository/ It would be a great addition!

– Sjoerd Smit – 2019-03-15T16:44:04.833