# Registering the user-defined function "multivariate normal cdf" for nonlinear modeling in JuMP

**URL:** <https://discourse.julialang.org/t/registering-the-user-defined-function-multivariate-normal-cdf-for-nonlinear-modeling-in-jump/79825>\
**Category:** Optimization (Mathematical)\
**Tags:** jump, nonlinear\
**Created:** [April 22, 2022, 9:32am UTC](https://discourse.julialang.org/t/registering-the-user-defined-function-multivariate-normal-cdf-for-nonlinear-modeling-in-jump/79825 "2022-04-22T09:32:36Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![naouess](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/naouess/32/35619_2.png) [@naouess](https://discourse.julialang.org/u/naouess)\
**Post date:** [April 22, 2022, 9:32am UTC](https://discourse.julialang.org/t/registering-the-user-defined-function-multivariate-normal-cdf-for-nonlinear-modeling-in-jump/79825/1 "2022-04-22T09:32:36Z")

</div>

Hello Julia-community,

I would like to use the multivariate normal cumulative distribution function within a constraint in JuMP.

For the function, I’m using the `mvnormcdf` from this package [MvNormalCDF.jl](https://github.com/PharmCat/MvNormalCDF.jl). Since this requires nonlinear modeling, I have to define the gradient and register the function as per [the documentation here](https://jump.dev/JuMP.jl/stable/manual/nlp/#Multivariate-functions).

I have the mathematical formulations on how to define the gradient but I’m having issues to write this in Julia/JuMP. The point where I’m stuck now, is that the 'mvnormcdf` returns two values (probability and error), for my nonlinear constraint however I only need the probability value.

Here is a simpflied version of my code:

```julia
using JuMP, LinearAlgebra
using Distributions, NLopt
using Random
using MvNormalCDF

begin
    ex2 = Model(NLopt.Optimizer)
    set_optimizer_attribute(ex2, "algorithm", :LD_MMA)

    @variable(ex2, x[1:2] >= 0)

    @objective(ex2, Min, x[1] + x[2])

    set_start_value(x[1], 5)
    set_start_value(x[2], 5)

    p = 0.9 # probability

    μ = [5; 5] # mean
    vari = [4; 4] # variance
    ρ = 0.5 # correlation
    Σ = [vari[1] ρ; ρ vari[2]] # covariance matrix
    a = [-Inf; -Inf] 

    norm_cdf(x) = cdf(Normal(), x)
    norm_pdf(x) = exp(-(x^2) / 2) / sqrt(2 * pi)
    register(ex2, :cdf, 1, norm_cdf, norm_pdf; autodiff = true)

    # Define the mvnormcdf and its gradient, then register the function 
    prob(x, y) = mvnormcdf(Σ, a-μ, [x; y]-μ) 
    function ∇prob(g::AbstractVector{T}, x::T, y::T) where {T}
        g[1] = norm_pdf((x - μ[1])/sqrt(vari[1])) * cdf((y-(μ[2]+ρ/sqrt(vari[1])*(x-μ[1])))/sqrt(vari[2]^2*(1- ρ^2)))
        g[2] = norm_pdf((y - μ[2])/sqrt(vari[2]))* cdf((x-(μ[1]+ρ/sqrt(vari[2])*(y-μ[2])))/sqrt(vari[1]^2*(1- ρ^2)))
        return
    end
    register(ex2, :mvncdf, 2, prob, ∇prob)

    # Use the mvncdf in the nonlinear constraint
    # This is where the type error occurs
    @NLconstraint(ex2, mvncdf(x[1], x[2]) >= p)

    JuMP.optimize!(ex2)

    println(value.(x))

end

```

The error I get is that the optimizer is getting a tuple instead of a float

```julia
TypeError: in typeassert, expected Float64, got a value of type Tuple{Float64,Float64}
in top-level scope at example.jl:42
in optimize! at JuMP/e0Uc2/src/optimizer_interface.jl:106
in optimize! at JuMP/e0Uc2/src/optimizer_interface.jl:106 
in #optimize!#100 at JuMP/e0Uc2/src/optimizer_interface.jl:130
in optimize! at MathOptInterface/YDdD3/src/Utilities/cachingoptimizer.jl:252
in optimize! at MathOptInterface/YDdD3/src/Bridges/bridge_optimizer.jl:319
in optimize! at NLopt/9omMe/src/MOI_wrapper.jl:942
in at NLopt/9omMe/src/MOI_wrapper.jl:907
in eval_constraint_jacobian at JuMP/e0Uc2/src/nlp.jl:583
in macro expansion at base/timing.jl:233 
in macro expansion at JuMP/e0Uc2/src/nlp.jl:585 
in _forward_eval_all at JuMP/e0Uc2/src/nlp.jl:503
in forward_eval at JuMP/e0Uc2/src/_Derivatives/forward.jl:163
in eval_and_check_return_type at JuMP/e0Uc2/src/_Derivatives/forward.jl:5
in eval_objective at JuMP/e0Uc2/src/nlp.jl:1169

```

Any hints on how to solve this?

---

<div class="post-metadata">

**Author:** ![jd-foster](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jd-foster/32/35824_2.png) [@jd-foster](https://discourse.julialang.org/u/jd-foster)\
**Post date:** [April 22, 2022, 11:06am UTC](https://discourse.julialang.org/t/registering-the-user-defined-function-multivariate-normal-cdf-for-nonlinear-modeling-in-jump/79825/2 "2022-04-22T11:06:20Z")

</div>

Welcome. Moving this over to the Optimization topic area.

---

<div class="post-metadata">

**Author:** ![blob](https://avatars.discourse-cdn.com/v4/letter/b/ebca7d/32.png) [@blob](https://discourse.julialang.org/u/blob)\
**Post date:** [April 22, 2022, 11:22am UTC](https://discourse.julialang.org/t/registering-the-user-defined-function-multivariate-normal-cdf-for-nonlinear-modeling-in-jump/79825/3 "2022-04-22T11:22:26Z")

</div>

You can set up `prob(x,y)` as the first value returned from `mvnormcdf`:

```julia
prob(x, y) = mvnormcdf(Σ, a-μ, [x; y]-μ)[1]

```

Also, your `cdf` in `∇prob` need the first argument to work, e.g.:

```julia
cdf(Normal(),(y-(μ[2]+ρ/sqrt(vari[1])*(x-μ[1])))/sqrt(vari[2]^2*(1- ρ^2))) #I added Normal(), not sure if this is correct for you

```

Then this works:

```julia

using JuMP, LinearAlgebra
using Distributions, NLopt
using Random
using MvNormalCDF

begin
    ex2 = Model(NLopt.Optimizer)
    set_optimizer_attribute(ex2, "algorithm", :LD_MMA)

    @variable(ex2, x[1:2] >= 0)

    @objective(ex2, Min, x[1] + x[2])

    set_start_value(x[1], 5)
    set_start_value(x[2], 5)

    p = 0.9 # probability

    μ = [5; 5] # mean
    vari = [4; 4] # variance
    ρ = 0.5 # correlation
    Σ = [vari[1] ρ; ρ vari[2]] # covariance matrix
    a = [-Inf; -Inf] 

    norm_cdf(x) = cdf(Normal(), x)
    norm_pdf(x) = exp(-(x^2) / 2) / sqrt(2 * pi)
    register(ex2, :cdf, 1, norm_cdf, norm_pdf; autodiff = true)

    # Define the mvnormcdf and its gradient, then register the function 
    prob(x, y) = mvnormcdf(Σ, a-μ, [x; y]-μ)[1]
    function ∇prob(g::AbstractVector{T}, x::T, y::T) where {T}
        g[1] = norm_pdf((x - μ[1])/sqrt(vari[1])) * cdf(Normal(),(y-(μ[2]+ρ/sqrt(vari[1])*(x-μ[1])))/sqrt(vari[2]^2*(1- ρ^2)))
        g[2] = norm_pdf((y - μ[2])/sqrt(vari[2]))* cdf(Normal(),(x-(μ[1]+ρ/sqrt(vari[2])*(y-μ[2])))/sqrt(vari[1]^2*(1- ρ^2)))
        return
    end

    register(ex2, :mvncdf, 2, prob, ∇prob)

    # Use the mvncdf in the nonlinear constraint
    # This is where the type error occurs
    @NLconstraint(ex2,mvncdf(x[1], x[2]) >= p)

    JuMP.optimize!(ex2)

    println(value.(x))

end

```

And I got:

```julia
[8.247724657052249, 8.247724657052249]

```

---

<div class="post-metadata">

**Author:** ![naouess](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/naouess/32/35619_2.png) [@naouess](https://discourse.julialang.org/u/naouess)\
**Post date:** [April 22, 2022, 11:46am UTC](https://discourse.julialang.org/t/registering-the-user-defined-function-multivariate-normal-cdf-for-nonlinear-modeling-in-jump/79825/4 "2022-04-22T11:46:16Z")

</div>

Thank you! This worked.