From 9b2b0048c5b1955d16d42b507b8c7f34f757e569 Mon Sep 17 00:00:00 2001
From: Andreas Noack <andreas@noack.dk>
Date: Thu, 10 Nov 2022 11:38:07 +0100
Subject: [PATCH 1/2] Test source with JET

This revealed some bugs in the untested bivariate T cdf. I've fixed
the bugs and also added a few tests.
---
 Project.toml     |  3 ++-
 src/tvpack.jl    | 22 ++++++++++++----------
 test/runtests.jl |  5 +++++
 test/tvpack.jl   | 11 +++++++++++
 4 files changed, 30 insertions(+), 11 deletions(-)

diff --git a/Project.toml b/Project.toml
index a69a96d..50f0009 100644
--- a/Project.toml
+++ b/Project.toml
@@ -26,8 +26,9 @@ julia = "1.3"
 [extras]
 ChainRulesTestUtils = "cdddcdb0-9152-4a09-a978-84456f9df70a"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
+JET = "c3a54625-cd67-489e-a8e7-0a5a0ff4e31b"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["ChainRulesTestUtils", "ForwardDiff", "Random", "Test"]
+test = ["ChainRulesTestUtils", "ForwardDiff", "JET", "Random", "Test"]
diff --git a/src/tvpack.jl b/src/tvpack.jl
index 7fdcd15..d20ac1e 100644
--- a/src/tvpack.jl
+++ b/src/tvpack.jl
@@ -12,7 +12,7 @@
 # Original source available from
 # http://www.math.wsu.edu/faculty/genz/software/fort77/tvpack.f
 
-function tvtcdf(nu::Int, h::Vector{Float64}, r::Vector{Float64})
+function tvtcdf(nu::Int, h::NTuple{3,Float64}, r::NTuple{3,Float64})
     pt = 0.5 * pi
     h1 = h[1]
     h2 = h[2]
@@ -128,10 +128,10 @@ end
 # # One Dimensional Globally Adaptive Integration Function
 
 function adonet(f::Function, a::Float64, b::Float64, tol::Float64)
-    ai = Vector{Float64}(100)
-    bi = Vector{Float64}(100)
-    ei = Vector{Float64}(100)
-    fi = Vector{Float64}(100)
+    ai = Vector{Float64}(undef, 100)
+    bi = Vector{Float64}(undef, 100)
+    ei = Vector{Float64}(undef, 100)
+    fi = Vector{Float64}(undef, 100)
     ai[1] = a
     bi[1] = b
     err = 1.0
@@ -242,7 +242,9 @@ end
 
 function bvtcdf(nu::Int, dh::Float64, dk::Float64, r::Float64)
     if nu < 1
-        return bvnuppercdf(-dh, -dk, r)
+        # The source call the commented out line below. Instead we error out
+        throw(DomainError(nu, "degrees of freedom parameter must be positive"))
+        # return bvnuppercdf(-dh, -dk, r)
     elseif 1.0 - r < eps()
         return tcdf(nu, min(dh, dk))
     elseif r + 1.0 < eps()
@@ -266,12 +268,12 @@ function bvtcdf(nu::Int, dh::Float64, dk::Float64, r::Float64)
         hs = sign(hrk)
         ks = sign(krh)
         if mod(nu, 2) == 0
-            bvt = atan2(sqrt(ors), -r) / (2.0pi)
+            bvt = atan(sqrt(ors), -r) / (2.0pi)
             gmph = dh / sqrt(16.0 * (nu + dh^2))
             gmpk = dk / sqrt(16.0 * (nu + dk^2))
-            btnckh = 2.0atan2(sqrt(xnkh), sqrt(1.0 - xnkh)) / pi
+            btnckh = 2.0atan(sqrt(xnkh), sqrt(1.0 - xnkh)) / pi
             btpdkh = 2.0sqrt(xnkh * (1.0 - xnkh)) / pi
-            btnchk = 2.0atan2(sqrt(xnhk), sqrt(1.0 - xnhk)) / pi
+            btnchk = 2.0atan(sqrt(xnhk), sqrt(1.0 - xnhk)) / pi
             btpdhk = 2.0sqrt(xnhk * (1.0 - xnhk)) / pi
             for j = 1:div(nu, 2)
                 bvt += gmph * (1.0 + ks * btnckh)
@@ -288,7 +290,7 @@ function bvtcdf(nu::Int, dh::Float64, dk::Float64, r::Float64)
             hkrn = dh*dk + r*nu
             hkn = dh*dk - nu
             hpk = dh + dk
-            bvt = atan2(-snu * (hkn * qhrk + hpk * hkrn), hkn * hkrn - nu * hpk * qhrk) / (2.0pi)
+            bvt = atan(-snu * (hkn * qhrk + hpk * hkrn), hkn * hkrn - nu * hpk * qhrk) / (2.0pi)
             if bvt < -eps()
                 bvt += 1.0
             end
diff --git a/test/runtests.jl b/test/runtests.jl
index 9534e00..26e67e6 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -5,3 +5,8 @@ for t in tests
     println("* running $fp")
     include(fp)
 end
+
+using JET
+@testset "Static analysis" begin
+    @test isempty(JET.get_reports(report_package("StatsFuns", target_modules=(StatsFuns,))))
+end
diff --git a/test/tvpack.jl b/test/tvpack.jl
index 66adbd1..fa4a6e0 100644
--- a/test/tvpack.jl
+++ b/test/tvpack.jl
@@ -1,4 +1,5 @@
 using StatsFuns, Test
+using StatsFuns: bvtcdf
 
 @testset "bvncdf: bivariate normal cdf" begin
 
@@ -28,3 +29,13 @@ using StatsFuns, Test
         @test !isnan(bvncdf(x,y,r))
     end
 end
+
+@testset "bvtcdf: bivariate T CDF" begin
+    # Tested against R's mvtnorm::pmvt
+    @test_throws DomainError bvtcdf(-1, 1.0, 1.0, 0.5)
+    @test bvtcdf(2, 1.0, 1.0, 1.0) ≈ 0.7886751345948129
+    @test bvtcdf(2, 1.0, 1.0, -1.0) ≈ 0.5773502691896257
+    @test bvtcdf(2, 1.0, -1.0, -1.0) == 0.0
+    @test bvtcdf(2, 1.0, 1.0, 0.5) ≈ 0.6811385938470963
+    @test bvtcdf(3, 1.0, 1.0, 0.5) ≈ 0.7001771920356403
+end

From 472871651d1aa0650f5600c4df2ae424affec2ae Mon Sep 17 00:00:00 2001
From: Andreas Noack <andreas@noack.dk>
Date: Thu, 10 Nov 2022 12:57:23 +0100
Subject: [PATCH 2/2] Update test/runtests.jl

Co-authored-by: David Widmann <devmotion@users.noreply.github.com>
---
 test/runtests.jl | 8 +++++---
 1 file changed, 5 insertions(+), 3 deletions(-)

diff --git a/test/runtests.jl b/test/runtests.jl
index 26e67e6..be71501 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -6,7 +6,9 @@ for t in tests
     include(fp)
 end
 
-using JET
-@testset "Static analysis" begin
-    @test isempty(JET.get_reports(report_package("StatsFuns", target_modules=(StatsFuns,))))
+if VERSION >= v"1.7"
+    using JET
+    @testset "Static analysis" begin
+        @test isempty(JET.get_reports(report_package("StatsFuns", target_modules=(StatsFuns,))))
+    end
 end