# JuMP.jl user-defined function has to return a scalar?

**URL:** <https://discourse.julialang.org/t/jump-jl-user-defined-function-has-to-return-a-scalar/11202>\
**Category:** Optimization (Mathematical)\
**Tags:** question, jump\
**Created:** [May 28, 2018, 3:11pm UTC](https://discourse.julialang.org/t/jump-jl-user-defined-function-has-to-return-a-scalar/11202 "2018-05-28T15:11:06Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![floswald](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/floswald/32/195_2.png) [@floswald](https://discourse.julialang.org/u/floswald)\
**Post date:** [May 28, 2018, 3:11pm UTC](https://discourse.julialang.org/t/jump-jl-user-defined-function-has-to-return-a-scalar/11202/1 "2018-05-28T15:11:06Z")

</div>

I have a question related to what was asked [here](https://discourse.julialang.org/t/interpolation-in-constraint/7325) some time ago. I want to use JuMP and need to differentiate through a function from the `BasisMatrices.jl` package and I am having trouble making that work. I don’t understand why the docs of JuMP say that the package uses ForwardDiff.jl (my error contains calls to ReverseSparseDiff instead). The call to `ForwardDiff.jl` works. artificial example:

```julia
module trying

	using JuMP
	using Ipopt
	using BasisMatrices
	using ForwardDiff

	function t1()
		chpar = ChebParams(10, 0.2, 2.0)
		Ψprime(k) = evalbase(chpar, k); # evalbase returns (1,10) matrix
		f(x) = ForwardDiff.jacobian(Ψprime, x);
		f([1.0])
	end

	function t2()
		chpar = ChebParams(10, 0.2, 2.0)
		m = Model(solver=IpoptSolver())
		@variable(m,x[1:3])

		function myfun(x)# evalbase returns (1,10) matrix
			return BasisMatrices.evalbase(chpar,x)
		end
		
		JuMP.register(m,:myfun,1,myfun,autodiff=true)

		@NLconstraint(m,sum(myfun(x[i]) for i=1:3) <= 0)
		@NLobjective(m,Min,1)
		solve(m)
	end
end

julia> trying.t1()
10×1 Array{Float64,2}:
  0.0     
  1.11111 
 -0.493827
 -3.16872 
  1.92654 
  4.74606 
 -4.15542 
 -5.56715 
  6.95444 
  5.41916 

julia> trying.t2()
ERROR: TypeError: #forward_eval#7: in typeassert, expected Float64, got Array{Float64,2}
Stacktrace:
 [1] #forward_eval#7(::ReverseDiffSparse.UserOperatorRegistry, ::Function, ::Array{Float64,1}, ::Array{Float64,1}, ::Array{ReverseDiffSparse.NodeData,1}, ::SparseMatrixCSC{Bool,Int64}, ::Array{Float64,1}, ::Array{Float64,1}, ::Array{Float64,1}, ::Array{Float64,1}, ::Array{Float64,1}, ::Array{Float64,1}) at /Users/74097/.julia/v0.6/ReverseDiffSparse/src/forward.jl:149
 [2] (::ReverseDiffSparse.#kw##forward_eval)(::Array{Any,1}, ::ReverseDiffSparse.#forward_eval, ::Array{Float64,1}, ::Array{Float64,1}, ::Array{ReverseDiffSparse.NodeData,1}, ::SparseMatrixCSC{Bool,Int64}, ::Array{Float64,1}, ::Array{Float64,1}, ::Array{Float64,1}, ::Array{Float64,1}, ::Array{Float64,1}, ::Array{Float64,1}) at ./<missing>:0
 [3] forward_eval_all(::JuMP.NLPEvaluator, ::Array{Float64,1}) at /Users/74097/.julia/v0.6/JuMP/src/nlp.jl:450
 [4] eval_hesslag(::JuMP.NLPEvaluator, ::Array{Float64,1}, ::Array{Float64,1}, ::Float64, ::Array{Float64,1}) at /Users/74097/.julia/v0.6/JuMP/src/nlp.jl:754
 [5] initialize(::JuMP.NLPEvaluator, ::Array{Symbol,1}) at /Users/74097/.julia/v0.6/JuMP/src/nlp.jl:399
 [6] loadproblem!(::Ipopt.IpoptMathProgModel, ::Int64, ::Int64, ::Array{Float64,1}, ::Array{Float64,1}, ::Array{Float64,1}, ::Array{Float64,1}, ::Symbol, ::JuMP.NLPEvaluator) at /Users/74097/.julia/v0.6/Ipopt/src/IpoptSolverInterface.jl:44
 [7] _buildInternalModel_nlp(::JuMP.Model, ::JuMP.ProblemTraits) at /Users/74097/.julia/v0.6/JuMP/src/nlp.jl:1248
 [8] #build#119(::Bool, ::Bool, ::JuMP.ProblemTraits, ::Function, ::JuMP.Model) at /Users/74097/.julia/v0.6/JuMP/src/solvers.jl:300
 [9] (::JuMP.#kw##build)(::Array{Any,1}, ::JuMP.#build, ::JuMP.Model) at ./<missing>:0
 [10] #solve#116(::Bool, ::Bool, ::Bool, ::Array{Any,1}, ::Function, ::JuMP.Model) at /Users/74097/.julia/v0.6/JuMP/src/solvers.jl:168
 [11] t2() at /Users/74097/git/HedonicRentals/julia/test2.jl:30

```

thanks for any help!

---

<div class="post-metadata">

**Author:** ![Juser](https://avatars.discourse-cdn.com/v4/letter/j/34f0e0/32.png) [@Juser](https://discourse.julialang.org/u/Juser)\
**Post date:** [May 29, 2018, 2:41am UTC](https://discourse.julialang.org/t/jump-jl-user-defined-function-has-to-return-a-scalar/11202/2 "2018-05-29T02:41:19Z")

</div>

I believe that is correct.

If I understand correctly, you are trying to create 10 constraints (that the element-wise sum of myfun(x[i]) for i=1:3 is elementwise non-positive. Such vectorized constraints only work with linear and quadratic constraints ([Link](http://www.juliaopt.org/JuMP.jl/0.18/refexpr.html#vectorized-operations)).

> Finally, note that this feature is not currently supported directly in nonlinear expressions; for example, a matrix–vector product will not work inside a call to the @NLconstraint macro.

You can still create such constraints element-wise, but it will involve making 10x more call to myfun than you were hoping.

---

<div class="post-metadata">

**Author:** ![MaximilianJHuber](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/maximilianjhuber/32/2579_2.png) [@MaximilianJHuber](https://discourse.julialang.org/u/MaximilianJHuber)\
**Post date:** [May 31, 2018, 3:48pm UTC](https://discourse.julialang.org/t/jump-jl-user-defined-function-has-to-return-a-scalar/11202/3 "2018-05-31T15:48:15Z")

</div>

My workaround back then was to code up a function that evaluates the `n`th-order Chebyshev at `x` in the domain:

```julia
function my_cheb(x, n)
    k_stst = 4.628988089138438
    a = 0.2*k_stst
    b = 2*k_stst
        
    z = (2/(b-a)) * (x-(a+b)/2)
    
    if (n == 1)
        return 1.
    elseif (n == 2)
        return z
    else
        
        out1 = 1.
        out2 = z
        out3 = 0.

        for j in 3:n
            out3 = 2. * z * out2 - out1
            out1 = out2
            out2 = out3
        end

        return out3
    end
end

```

This, of course, involves repeated computation of already known numbers. And I am uncertain about type stability when JuMP calls the function.

It is a bit better to hand the Chebyshev coefficients to the user-defined function, like:

```julia
function Ψprime(k, θ...)
    (BasisMatrix(Basis(ChebParams(10, 0.2*4.628988089138438, 2*4.628988089138438)), 
        Expanded(), [k]).vals[1] * collect(θ))[1]
end

```

It returns a scalar too, but it saves some computations compared to the solution above. The function call inside some `@NLconstraint` is then:

```julia
Ψprime(Kprime[i], θ[1], θ[2], θ[3], θ[4], θ[5], θ[6], θ[7], θ[8], θ[9], θ[10])

```