13
$\begingroup$

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.

$\endgroup$
  • 3
    $\begingroup$ They compile for real values, but do not appear to compile for complex arguments. $\endgroup$ – asim Feb 20 '13 at 17:56
  • 1
    $\begingroup$ 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. $\endgroup$ – Oleksandr R. Feb 20 '13 at 18:25
  • 1
    $\begingroup$ On the other hand, regarding the PDF of the Voigt distribution: dx.doi.org/10.1016/0368-2048(94)02189-7 $\endgroup$ – Oleksandr R. Feb 20 '13 at 18:46
7
$\begingroup$

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]}]
       ]
$\endgroup$
9
$\begingroup$

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.

$\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.