Merge pull request #31 from JoelTrent/dev

Increment version to 0.1.2.

Project.toml

@@ -1,12 +1,11 @@
 name = "EllipseSampling"
 uuid = "7d220c50-6298-48cd-9710-1d1a8ef13806"
 authors = ["JoelTrent <[email protected]> and contributors"]
-version = "0.1.1"
+version = "0.1.2"
 
 [deps]
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 Elliptic = "b305315f-e792-5b7a-8f41-49f472929428"
-LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665"
 SnoopPrecompile = "66db9d55-30c0-4569-8b51-7e840670fc0c"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
@@ -20,7 +19,8 @@ StaticArrays = "1.3"
 julia = "1.1"
 
 [extras]
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Test"]
+test = ["LinearAlgebra", "Test"]

src/matrix_to_parameters.jl

@@ -23,11 +23,26 @@ function calculate_ellipse_parameters(Γ::Matrix{Float64}, ind1::Int, ind2::Int,
     (0 < ind1 && ind1 ≤ size(Γ)[1]) || throw(BoundsError("ind1 must be a valid row index in Γ.")) 
     (0 < ind2 && ind2 ≤ size(Γ)[1]) || throw(BoundsError("ind2 must be a valid row index in Γ."))
 
-    Hw = inv(Γ[[ind1, ind2], [ind1, ind2]]) .* 0.5 ./ (Distributions.quantile(Distributions.Chisq(2), confidence_level)*0.5) # normalise Hw so that the RHS of the ellipse equation == 1
+    # normalise Hw so that the RHS of the ellipse equation == 1
+    Hw = inv(Γ[[ind1, ind2], [ind1, ind2]]) .* 0.5 ./ (Distributions.quantile(Distributions.Chisq(2), confidence_level)*0.5)
 
     α = atan(2*Hw[1,2]/(Hw[1,1]-Hw[2,2]))/2
-    y_radius = sqrt( (cos(α)^4 - sin(α)^4) / (Hw[2,2]*(cos(α)^2) - Hw[1,1]*(sin(α)^2))  )
-    x_radius = sqrt( (cos(α)^2) / (Hw[1,1] - (sin(α)^2)/y_radius^2))
+
+    # if α close to +/-0.25pi, +/-1.25pi, then switch to BigFloat precision
+    if isapprox(abs(rem(α/pi, 1)), 0.25, atol=1e-2)
+
+        # convert values to BigFloat for enhanced precision - required for correct results when α → 0.25pi or 1.25pi.
+        Hw = inv(BigFloat.(Γ[[ind1, ind2], [ind1, ind2]], RoundUp, precision=64)) .* 0.5 ./ (Distributions.quantile(Distributions.Chisq(2), confidence_level)*0.5) 
+
+        α = atan(2*Hw[1,2]/(Hw[1,1]-Hw[2,2]))/2
+        y_radius = sqrt( (cos(α)^4 - sin(α)^4) / (Hw[2,2]*(cos(α)^2) - Hw[1,1]*(sin(α)^2))  )
+        x_radius = sqrt( (cos(α)^2) / (Hw[1,1] - (sin(α)^2)/y_radius^2))
+
+        α, x_radius, y_radius = convert(Float64, α), convert(Float64, x_radius), convert(Float64, y_radius)
+    else
+        y_radius = sqrt( (cos(α)^4 - sin(α)^4) / (Hw[2,2]*(cos(α)^2) - Hw[1,1]*(sin(α)^2))  )
+        x_radius = sqrt( (cos(α)^2) / (Hw[1,1] - (sin(α)^2)/y_radius^2))
+    end
 
     a = max(x_radius, y_radius)
     b = min(x_radius, y_radius)
@@ -41,6 +56,5 @@ function calculate_ellipse_parameters(Γ::Matrix{Float64}, ind1::Int, ind2::Int,
     # atan(eigs[2,1], eigs[1,1]) # if a is x axis
     # atan(eigs[2,1], eigs[1,1]) + pi/2 # if a is y axis
 
-
     return a, b, x_radius, y_radius, α 
 end

test/runtests.jl

@@ -137,5 +137,24 @@ end
         else
             @test isapprox_ellipsesampling(atan(eigvectors[2,1], eigvectors[1,1]) + 0.5*pi, α)
         end
+
+        # For issue #30 when α → +/- 0.25 or +/- 1.25
+        a, b = 2.0, 1.0 
+        α = 0.25*π
+
+        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))
+        Γ = convert.(Float64, inv(BigFloat.(Hw, precision=64)))
+
+        a_calc, b_calc, _, _, _ = EllipseSampling.calculate_ellipse_parameters(Γ, 1, 2, confidence_level)
+
+        @test isapprox_ellipsesampling(a, a_calc)
+        @test isapprox_ellipsesampling(b, b_calc)
     end
+
 end