`calculate_ellipse_parameters` fails to calculate the correct values for `a` and `b` when the ellipse rotation approaches 0.25 or 1.25

[JoelTrent]

Minimum working example (for version 0.1.1):

using EllipseSampling
import Distributions
using Plots
gr()

N=8

a, b = 2.0, 1.0 
α = 0.25*π
Cx, Cy= 2.0, 2.0
points=generate_N_equally_spaced_points(N, a, b, α, Cx, Cy, start_point_shift=0.0)
plt=scatter(points[1,:], points[2,:], label=nothing, mw=0, ms=8, aspect_ratio=:equal)

Hw11 = (cos(α)^2 / a^2 + sin(α)^2 / b^2)
Hw22 = (sin(α)^2 / a^2 + cos(α)^2 / b^2)
Hw12 = cos(α)*sin(α)*(1/a^2 - 1/b^2)
Hw_norm = [Hw11 Hw12; Hw12 Hw22]

confidence_level=0.95
Hw = Hw_norm ./ (0.5 ./ (Distributions.quantile(Distributions.Chisq(2), confidence_level)*0.5))
Γ = inv(Hw)

points2 = generate_N_equally_spaced_points(N, Γ, [Cx, Cy], 1, 2, confidence_level=confidence_level, start_point_shift=0.0)
scatter!(points2[1, :], points2[2, :], label=nothing, markeralpha=0.5, mw=0, ms=4)

display(plt)

If we modify the rotation slightly away from 0.25pi the problem goes away:

a, b = 2.0, 1.0 
α = 0.251*π
Cx, Cy= 2.0, 2.0
points=generate_N_equally_spaced_points(N, a, b, α, Cx, Cy, start_point_shift=0.0)
plt=scatter(points[1,:], points[2,:], label=nothing, mw=0, ms=8, aspect_ratio=:equal)

Hw11 = (cos(α)^2 / a^2 + sin(α)^2 / b^2)
Hw22 = (sin(α)^2 / a^2 + cos(α)^2 / b^2)
Hw12 = cos(α)*sin(α)*(1/a^2 - 1/b^2)
Hw_norm = [Hw11 Hw12; Hw12 Hw22]

confidence_level=0.95
Hw = Hw_norm ./ (0.5 ./ (Distributions.quantile(Distributions.Chisq(2), confidence_level)*0.5))
Γ = inv(Hw)

points2 = generate_N_equally_spaced_points(N, Γ, [Cx, Cy], 1, 2, confidence_level=confidence_level, start_point_shift=0.0)
scatter!(points2[1, :], points2[2, :], label=nothing, markeralpha=0.5, mw=0, ms=4)

display(plt)

We can confirm this issue by inspecting the ellipse values returned by calculate_ellipse_parameters:

# for α = 0.25*π
a_calc, b_calc, _, _, α_calc = EllipseSampling.calculate_ellipse_parameters(Γ, 1, 2, confidence_level)

println("a=", a_calc); println("b=", b_calc); println("α/pi=", α_calc/pi)

a=1.7994708216848752
b=1.0307764064044151
α=0.24999999999999992

[JoelTrent]

The culprit is likely to be line 29 in src/matrix_to_parameters.jl which subtracts cos(α)^4 - sin(α)^4. When alpha is 0.25pi and 1.25pi cos(α) is equal to sin(α). However due to numerical issues, both when calculating alpha and when calculating the subtraction of the 4th power of these functions, the actual value is not quite a zero:

println(cos(0.25pi)^4 - sin(0.25pi)^4); println(cos(α)^4 - sin(α)^4)

1.3877787807814457e-16
1.3877787807814457e-16

The behaviour has also been noted for negative values of alpha close to -0.25pi and -1.25pi.

We can fix this issue by enhancing the precision of our calculations using BigFloat on the normalised matrix Hw on line 26. A precision of 64 gives a satisfactory result, however, it does impact the performance of the function by about 40 times. Lower precision values also gave good results for this test case, however, the additional precision was valuable when the magnitude of the two variables is significantly different.

A satisfactory compromise is to selectively enhance the precision if alpha is detected to be close to these problematic values using: isapprox(abs(rem(α/pi, 1)), 0.25, atol=1e-2)