Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

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.

share|improve this question
3  
They compile for real values, but do not appear to compile for complex arguments. – asim Feb 20 '13 at 17:56
1  
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
1  
On the other hand, regarding the PDF of the Voigt distribution: dx.doi.org/10.1016/0368-2048(94)02189-7 – Oleksandr R. Feb 20 '13 at 18:46
up vote 4 down vote accepted

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]}]
       ]
share|improve this answer

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.

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.