0.3 Less allocation, ForwardDiff compat, #8

Project.toml

@@ -3,7 +3,7 @@ uuid = "37188c8d-bc69-4638-b057-733e744175ec"
 keywords = ["cdf", "multivariate-normal", "distribution"]
 desc = "Numerical Computation of Multivariate Normal Probabilities."
 authors = ["PharmCat <v.s.arnautov@yandex.ru>", "Andrew Gough"]
-version = "0.2.6"
+version = "0.3.0"
 
 
 [deps]
@@ -23,8 +23,9 @@ julia = "1"
 
 [extras]
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
+ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
 
 [targets]
-test = ["Test", "Distributions", "StableRNGs"]
+test = ["Test", "Distributions", "ForwardDiff", "StableRNGs"]

README.md

@@ -86,7 +86,7 @@ m = 5000
 ```
 
 ```
-     mvnormcdf(Σ::AbstractMatrix{<:Real}, a::AbstractVector{<:Real}, b::AbstractVector{<:Real}; m::Integer = 1000*size(Σ,1), rng = RandomDevice())
+     mvnormcdf(Σ::AbstractMatrix, a::AbstractVector, b::AbstractVector; m::Integer = 1000*size(Σ,1), rng = RandomDevice())
 ```
 
 Computes the Multivariate Normal probability integral with covariance matrix Σ, mean [0,...].
@@ -101,6 +101,14 @@ b = [2; 2]
 #(0.4306346895870772, 0.00015776288569406053)
 ```
 
+## Not exported
+
+```
+    MvNormalCDF.qsimvnv!(Σ::AbstractMatrix{T}, a::AbstractVector{T}, b::AbstractVector{T}, m::Integer, rng) where T
+```
+
+Re-coded in Julia from the MATLAB function qsimvnv(m,r,a,b); mutate Σ, a, b. Computes the Multivariate Normal probability integral with covariance matrix Σ, mean [0,...].
+
 # Reference
 - Genz, A. (1992). Numerical computation of multivariate normal probabilities. Journal of Computational and Graphical Statistics, 1, 141--150
 - Genz, A. (1993). Comparison of methods for the computation of multivariate normal probabilities. Computing Science and Statistics, 25, 400--405

src/functions.jl

@@ -72,7 +72,7 @@ function mvnormcdf(dist::MvNormal, a, b; m::Integer = 1000*size(dist.Σ,1), rng
 end
 
 """
-    mvnormcdf(μ::AbstractVector{<:Real}, Σ::AbstractMatrix{<:Real}, a::AbstractVector{<:Real}, b::AbstractVector{<:Real}; m::Integer = 1000*size(Σ,1), rng = RandomDevice())
+    mvnormcdf(μ::AbstractVector, Σ::AbstractMatrix, a::AbstractVector, b::AbstractVector; m::Integer = 1000*size(Σ,1), rng = RandomDevice())
 
 Computes the Multivariate Normal probability integral using a quasi-Monte-Carlo
 algorithm with m points for positive definite covariance matrix Σ, mean [0,...], with lower
@@ -103,12 +103,12 @@ m = 5000
 Results will vary slightly from run-to-run due to the quasi-Monte-Carlo
 algorithm.
 """
-function mvnormcdf(Σ::AbstractMatrix{<:Real}, a::AbstractVector{<:Real}, b::AbstractVector{<:Real}; m::Integer = 1000*size(Σ,1), rng = RandomDevice())
-    mvnormcdf(Zeros(size(Σ, 1)), Σ, a, b, m = m, rng = rng)
+function mvnormcdf(Σ::AbstractMatrix, a::AbstractVector, b::AbstractVector; m::Integer = 1000*size(Σ,1), rng = RandomDevice())
+    mvnormcdf(Zeros{promote_type(eltype(Σ), eltype(a), eltype(b), Float64)}(size(Σ, 1)), Σ, a, b, m = m, rng = rng)
 end
 
 """
-     mvnormcdf(Σ::AbstractMatrix{<:Real}, a::AbstractVector{<:Real}, b::AbstractVector{<:Real}; m::Integer = 1000*size(Σ,1), rng = RandomDevice())
+     mvnormcdf(Σ::AbstractMatrix, a::AbstractVector, b::AbstractVector; m::Integer = 1000*size(Σ,1), rng = RandomDevice())
 
 Non-central MVN distributions (with non-zero mean) can use this function by adjusting
 the integration limits. Subtract the mean vector, μ, from each
@@ -123,7 +123,8 @@ b = [2; 2]
 #(0.4306346895870772, 0.00015776288569406053)
 ```
 """
-function mvnormcdf(μ, Σ::AbstractMatrix{<:Real}, a::AbstractVector{<:Real}, b::AbstractVector{<:Real}; m::Integer = 1000*size(Σ,1), rng = RandomDevice())
+function mvnormcdf(μ::AbstractVector, Σ::AbstractMatrix, a::AbstractVector, b::AbstractVector; m::Integer = 1000*size(Σ,1), rng = RandomDevice()) 
+    T = promote_type(eltype(μ), eltype(Σ), eltype(a), eltype(b), Float64)
     # check for proper dimensions
     n  = size(Σ, 1)
     nc = size(Σ, 2) 	# assume square Cov matrix nxn
@@ -133,26 +134,26 @@ function mvnormcdf(μ, Σ::AbstractMatrix{<:Real}, a::AbstractVector{<:Real}, b:
     # check dimensions of lower vector, upper vector, and cov matrix match
     (n == length(a) == length(b)) || throw(DimensionMismatch("iconsistent argument dimensions. Sizes: Σ $(size(Σ))  a $(size(a))  b $(size(b))"))
 
-    # check that lower integration limit a < upper integration limit b for all elements
-    all(a .<= b) || throw(ArgumentError("lower integration limit a must be <= upper integration limit b"))
+    # check that lower integration limit a <= upper integration limit b for all elements
+    all( x -> x[1] <= x[2], zip(a, b)) || throw(ArgumentError("lower integration limit a must be <= upper integration limit b"))
     # check that Σ is positive definate; if not, print warning
     # isposdef(Σ) || @warn "covariance matrix Σ fails positive definite check"
     # check if Σ, a, or b contains NaNs
-    if any(isnan.(Σ)) || any(isnan.(a)) || any(isnan.(b))
+    if any(isnan, Σ) || any(isnan, a) || any(isnan, b)
         p = NaN
         e = NaN
         return (p, e)
     end
     # check if a==b
     if a == b
-        p = 0.0
-        e = 0.0
+        p = zero(T)
+        e = zero(T)
         return (p, e)
     end
     # check if a = -Inf & b = +Inf
-    if all(a .== -Inf) && all(b .== Inf)
-        p = 1.0
-        e = 0.0
+    if all(x -> x == -Inf, a) && all(x -> x == Inf, b)
+        p = one(T)
+        e = zero(T)
         return (p, e)
     end
     ##################################################################
@@ -161,7 +162,7 @@ function mvnormcdf(μ, Σ::AbstractMatrix{<:Real}, a::AbstractVector{<:Real}, b:
     # 3-dimesional Σ have exact solutions. Integration range [0,∞]
     #
     ##################################################################
-    if all(a .== zero(eltype(a))) && all(b .== Inf) && n <= 3
+    if all(iszero, a) && all(x -> x == Inf, b) && n <= 3
         #Σstd = sqrt.(diag(Σ))
         Σstd  = Vector{Float64}(undef, n)
         @inbounds for i in 1:n
@@ -175,21 +176,25 @@ function mvnormcdf(μ, Σ::AbstractMatrix{<:Real}, a::AbstractVector{<:Real}, b:
             p = 1/8 + (asin(Rcorr[1, 2]) + asin(Rcorr[2, 3]) + asin(Rcorr[1, 3])) / (4π)
             e = eps()
         end
-        return (p,e)
+        return (p, e)
     end
     #
-    a = a .- μ
-    b = b .- μ
-    qsimvnv!(copy_oftype(Σ, Float64), a, b, m, rng)
+    at = copy_oftype(a, T)
+    bt = copy_oftype(b, T)
+    at .-= μ
+    bt .-= μ
+    qsimvnv!(copy_oftype(Σ, T), at, bt, m, rng)
 end
 
 """
+    MvNormalCDF.qsimvnv!(Σ::AbstractMatrix{T}, a::AbstractVector{T}, b::AbstractVector{T}, m::Integer, rng) where T
 Re-coded in Julia from the MATLAB function qsimvnv(m,r,a,b)
 Alan Genz is the author the MATLAB qsimvnv() function.
 
 ! Mutate Σ, a, b.
 """
-function qsimvnv!(Σ::AbstractMatrix, a::AbstractVector{<:Real}, b::AbstractVector{<:Real}, m::Integer, rng)
+function qsimvnv!(Σ::AbstractMatrix{T}, a::AbstractVector{T}, b::AbstractVector{T}, m::Integer, rng) where T
+    #T = promote_type(T1, T2)
     ##################################################################
     #
     # get lower cholesky matrix and (potentially) re-ordered integration vectors
@@ -213,25 +218,25 @@ function qsimvnv!(Σ::AbstractMatrix, a::AbstractVector{<:Real}, b::AbstractVect
         if ai < 9ct
             c1 = cdf(ZDIST, ai / ct)
         else
-            c1 = 1.0
+            c1 = one(T)
         end
     else
-        c1 = 0.0
+        c1 = zero(T)
     end
     # if bi is +infinity, explicity set d=0
     if bi > -9ct
         if bi < 9ct
             d1 = cdf(ZDIST, bi / ct)
         else
-            d1 = 1.0
+            d1 = one(T)
         end
     else
-        d1 = 0.0
+        d1 = zero(T)
     end
     n   = size(Σ, 1) 	# assume square Cov matrix nxn
     cxi = c1			# initial cxi; genz uses ci but it conflicts with Lin. Alg. ci variable
     dci = d1 - cxi		# initial dcxi
-    p   = 0.0			# probablity = 0
+    p   = zero(T)       # probablity = 0
     e   = 0.0			# error = 0
     # Richtmyer generators
     ps  = sqrt.(primes(Int(floor(5 * n * log(n + 1) / 4)))) # Richtmyer generators
@@ -242,18 +247,18 @@ function qsimvnv!(Σ::AbstractMatrix, a::AbstractVector{<:Real}, b::AbstractVect
     Jnv    = ones(1, nv)
     cfill  = fill(cxi, nv) 	            # evaulate at nv quasirandom points row vec
     dpfill = fill(dci, nv)
-    y      = zeros(nv, n - 1)			# n-1 (cols), nv (rows), preset to zero # change row-col for col-operation
+    y      = zeros(T, nv, n - 1)			# n-1 (cols), nv (rows), preset to zero # change row-col for col-operation
     #=
     Randomization loop for ns samples
     j is the number of samples to integrate over,
     but each with a vector nv in length
     i is the number of dimensions, or integrals to comptue
     =#
-    c  = zeros(length(cfill))
-    dc = zeros(length(dpfill))
-    pv = zeros(length(dpfill))
-    d  = zeros(length(Jnv))
-    tv = zeros(nv)
+    c  = zeros(T, length(cfill))
+    dc = zeros(T, length(dpfill))
+    pv = zeros(T, length(dpfill))
+    d  = zeros(T, length(Jnv))
+    tv = zeros(T, nv)
 
     for j in 1:ns					# loop for ns samples
         copyto!(c, cfill)
@@ -336,7 +341,8 @@ Permutated lower input vector:
 bp = [1, 2, 4]
 """
 #############
-function _chlrdr!(Σ::AbstractMatrix{T}, a::AbstractVector{T}, b::AbstractVector{T}) where T <: Real
+function _chlrdr!(Σ::AbstractMatrix{T}, a::AbstractVector{T}, b::AbstractVector{T}) where T
+
     # define constants
     ep = 1e-10 # singularity tolerance
     ϵ  = eps()
@@ -348,6 +354,7 @@ function _chlrdr!(Σ::AbstractMatrix{T}, a::AbstractVector{T}, b::AbstractVector
     am  = zero(T)
     bm  = zero(T)
     im  = zero(T)
+    #c   = copyto!(Matrix{T}(undef, n, n), Σ)
     c   = Σ
     ap  = a
     bp  = b

test/runtests.jl

@@ -1,7 +1,7 @@
 # MvNormalCDF
 # Copyright © 2019-2021 Vladimir Arnautov aka PharmCat <mail@pharmcat.net>, Andrew Gough
 using MvNormalCDF
-using Test, Distributions, StableRNGs
+using Test, Distributions, StableRNGs, ForwardDiff
 
 td = Array{Any}(undef,(14,9))
 #  1-cov mtx 2-a 3-b 4-m 5-p 6-ptol 7-e 8-etol
@@ -259,6 +259,26 @@ td[14,1] = [59.227 2.601 3.38 8.303 -0.334 11.029 10.908 0.739 4.703 7.075 8.049
 
  end
 
+@testset "ForwardDiff test" begin
+  μ = [1., 2., 3.] 
+  Σ = [1 0.25 0.2; 0.25 1 0.333333333; 0.2 0.333333333 1]
+  ag = [-1; -4; -2]
+  bg = [1; 4; 2]
+  gf(x) = MvNormalCDF.mvnormcdf(x, Σ, ag, bg)[1]
+  @test_nowarn ForwardDiff.gradient(gf, [1, 2, 3])
+
+  μ = [1., 2., 3.] 
+  Σ = [1 0.25 0.2; 0.25 1 0.333333333; 0.2 0.333333333 1]
+  ag = [-1; -4; -2]
+  bg = [1; 4; 2]
+  gf2(x) = MvNormalCDF.mvnormcdf(μ, reshape(x, 3, 3), ag, bg)[1]
+  @test_nowarn ForwardDiff.gradient(gf2, [1, 0.25, 0.2, 0.25, 1, 0.333333333, 0.2, 0.333333333, 1])
+
+  f(x) = MvNormalCDF.mvnormcdf([0.;0.], reshape([0.1; 0.; 0.; 0.1], 2, 2), [x; -1.], [1.; 1.])[1]; 
+  ForwardDiff.derivative(f, -1.)
+end
+
+
 #=
 using BenchmarkTools
 mvn = MvNormal(td[14,1])
@@ -331,4 +351,4 @@ BenchmarkTools.Trial: 1525 samples with 1 evaluation.
 
  Memory estimate: 52.70 KiB, allocs estimate: 63.
 
-=#
+=#