# Automatic differentiation and Trapz

**URL:** <https://discourse.julialang.org/t/automatic-differentiation-and-trapz/82135>\
**Category:** Optimization (Mathematical)\
**Tags:** ad\
**Created:** [June 2, 2022, 2:06pm UTC](https://discourse.julialang.org/t/automatic-differentiation-and-trapz/82135 "2022-06-02T14:06:52Z")\
**Posts on this page:** 11\
**Page:** 1

<div class="post-metadata">

**Author:** ![chupin](https://avatars.discourse-cdn.com/v4/letter/c/f07891/32.png) [@chupin](https://discourse.julialang.org/u/chupin)\
**Post date:** [June 2, 2022, 2:06pm UTC](https://discourse.julialang.org/t/automatic-differentiation-and-trapz/82135/1 "2022-06-02T14:06:52Z")

</div>

It seems that Nonconvex is not able to deal with the trapz function, but maybe, I’ve done something wrong. The code (that is a nonsense toy code) :

```julia
using Nonconvex
Nonconvex.@load Ipopt
using Trapz

function obj(hV,N,M)
    h=reshape(hV,(N,M))
    vx=range(0,1,length=N)
    vy=range(0,1,length=M)
    output = trapz((vx,vy),h)
    return tr(h)
end

function cstr(hV,N,M,b)
    h=reshape(hV,(N,M))
    return det(h)-b
end

N=2
M=2

model = Model(x->obj(x,N,M))
addvar!(model, vec(zeros(N,M)), 4.0*vec(ones(N,M)))
add_ineq_constraint!(model, x -> cstr(x, N, M, 2.0))

alg = IpoptAlg()
options = IpoptOptions(first_order = true, tol = 1e-7)
r = optimize(model, alg, vec(rand(N,M)), options = options)

```

---

<div class="post-metadata">

**Author:** ![chupin](https://avatars.discourse-cdn.com/v4/letter/c/f07891/32.png) [@chupin](https://discourse.julialang.org/u/chupin)\
**Post date:** [June 14, 2022, 7:59am UTC](https://discourse.julialang.org/t/automatic-differentiation-and-trapz/82135/2 "2022-06-14T07:59:48Z")

</div>

Any idea ?

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [June 14, 2022, 8:47am UTC](https://discourse.julialang.org/t/automatic-differentiation-and-trapz/82135/3 "2022-06-14T08:47:34Z")

</div>

Cc @mohamed82008

What’s the error?

---

<div class="post-metadata">

**Author:** ![chupin](https://avatars.discourse-cdn.com/v4/letter/c/f07891/32.png) [@chupin](https://discourse.julialang.org/u/chupin)\
**Post date:** [June 14, 2022, 9:15am UTC](https://discourse.julialang.org/t/automatic-differentiation-and-trapz/82135/4 "2022-06-14T09:15:44Z")

</div>

```julia
ERROR: LoadError: Mutating arrays is not supported -- called setindex!(::Array{Float64, 0}, _...)
Stacktrace:
  [1] error(s::String)
    @ Base ./error.jl:33
  [2] (::Zygote.var"#437#438"{Array{Float64, 0}})(#unused#::Nothing)
    @ Zygote ~/.julia/packages/Zygote/DkIUK/src/lib/array.jl:71
  [3] (::Zygote.var"#2337#back#439"{Zygote.var"#437#438"{Array{Float64, 0}}})(Δ::Nothing)
    @ Zygote ~/.julia/packages/ZygoteRules/AIbCs/src/adjoint.jl:67
  [4] macro expansion
    @ ~/.julia/packages/Trapz/UoNri/src/kernels.jl:42 [inlined]
  [5] macro expansion
    @ ./simdloop.jl:77 [inlined]
  [6] Pullback
    @ ~/.julia/packages/Trapz/UoNri/src/kernels.jl:41 [inlined]
  [7] macro expansion
    @ ./none:0 [inlined]
  [8] Pullback
    @ ./none:0 [inlined]
  [9] Pullback
    @ ~/.julia/packages/Trapz/UoNri/src/Trapz.jl:7 [inlined]
 [10] (::typeof(∂(trapz)))(Δ::Nothing)
    @ Zygote ~/.julia/packages/Zygote/DkIUK/src/compiler/interface2.jl:0
 [11] Pullback
    @ ~/Nextcloud/Dauphine/OptiForm/Julia/Nonconvex/testTrapz.jl:9 [inlined]
 [12] (::typeof(∂(obj)))(Δ::Float64)
    @ Zygote ~/.julia/packages/Zygote/DkIUK/src/compiler/interface2.jl:0
 [13] Pullback
    @ ~/Nextcloud/Dauphine/OptiForm/Julia/Nonconvex/testTrapz.jl:21 [inlined]
 [14] (::typeof(∂(#1)))(Δ::Float64)
    @ Zygote ~/.julia/packages/Zygote/DkIUK/src/compiler/interface2.jl:0
 [15] #208
    @ ~/.julia/packages/Zygote/DkIUK/src/lib/lib.jl:207 [inlined]
 [16] (::Zygote.var"#1750#back#210"{Zygote.var"#208#209"{Tuple{Tuple{Nothing}}, typeof(∂(#1))}})(Δ::Float64)
    @ Zygote ~/.julia/packages/ZygoteRules/AIbCs/src/adjoint.jl:67
 [17] Pullback
    @ ~/.julia/packages/NonconvexCore/uZAVo/src/functions/functions.jl:170 [inlined]
 [18] (::typeof(∂(#_#8)))(Δ::Float64)
    @ Zygote ~/.julia/packages/Zygote/DkIUK/src/compiler/interface2.jl:0
 [19] (::Zygote.var"#208#209"{Tuple{Tuple{Nothing, Nothing}, Tuple{Nothing}}, typeof(∂(#_#8))})(Δ::Float64)
    @ Zygote ~/.julia/packages/Zygote/DkIUK/src/lib/lib.jl:207
 [20] #1750#back
    @ ~/.julia/packages/ZygoteRules/AIbCs/src/adjoint.jl:67 [inlined]
 [21] Pullback
    @ ~/.julia/packages/NonconvexCore/uZAVo/src/functions/functions.jl:170 [inlined]
 [22] (::typeof(∂(λ)))(Δ::Float64)
    @ Zygote ~/.julia/packages/Zygote/DkIUK/src/compiler/interface2.jl:0
 [23] Pullback
    @ ~/.julia/packages/NonconvexCore/uZAVo/src/models/vec_model.jl:79 [inlined]
 [24] (::typeof(∂(λ)))(Δ::Float64)
    @ Zygote ~/.julia/packages/Zygote/DkIUK/src/compiler/interface2.jl:0
 [25] #208
    @ ~/.julia/packages/Zygote/DkIUK/src/lib/lib.jl:207 [inlined]
 [26] (::Zygote.var"#1750#back#210"{Zygote.var"#208#209"{Tuple{Tuple{Nothing}}, typeof(∂(λ))}})(Δ::Float64)
    @ Zygote ~/.julia/packages/ZygoteRules/AIbCs/src/adjoint.jl:67
 [27] Pullback
    @ ~/.julia/packages/NonconvexCore/uZAVo/src/functions/functions.jl:170 [inlined]
 [28] (::typeof(∂(#_#8)))(Δ::Float64)
    @ Zygote ~/.julia/packages/Zygote/DkIUK/src/compiler/interface2.jl:0
 [29] (::Zygote.var"#208#209"{Tuple{Tuple{Nothing, Nothing}, Tuple{Nothing}}, typeof(∂(#_#8))})(Δ::Float64)
    @ Zygote ~/.julia/packages/Zygote/DkIUK/src/lib/lib.jl:207
 [30] #1750#back
    @ ~/.julia/packages/ZygoteRules/AIbCs/src/adjoint.jl:67 [inlined]
 [31] Pullback
    @ ~/.julia/packages/NonconvexCore/uZAVo/src/functions/functions.jl:170 [inlined]
 [32] (::typeof(∂(λ)))(Δ::Float64)
    @ Zygote ~/.julia/packages/Zygote/DkIUK/src/compiler/interface2.jl:0
 [33] Pullback
    @ ~/.julia/packages/NonconvexCore/uZAVo/src/functions/counting_function.jl:9 [inlined]
 [34] (::typeof(∂(λ)))(Δ::Float64)
    @ Zygote ~/.julia/packages/Zygote/DkIUK/src/compiler/interface2.jl:0
 [35] (::Zygote.var"#52#53"{typeof(∂(λ))})(Δ::Float64)
    @ Zygote ~/.julia/packages/Zygote/DkIUK/src/compiler/interface.jl:41
 [36] gradient(f::Function, args::Vector{Float64})
    @ Zygote ~/.julia/packages/Zygote/DkIUK/src/compiler/interface.jl:76
 [37] (::NonconvexIpopt.var"#eval_grad_f#13"{NonconvexCore.CountingFunction{NonconvexCore.Objective{NonconvexCore.var"#133#140"{Model{Vector{Any}}, NonconvexCore.Unflatten{Vector{Float64}, typeof(identity)}}, Base.RefValue{Float64}, Set{Symbol}}}})(x::Vector{Float64}, grad_f::Vector{Float64})
    @ NonconvexIpopt ~/.julia/packages/NonconvexIpopt/jkxtN/src/ipopt.jl:135
 [38] _Eval_Grad_F_CB(n::Int32, x_ptr::Ptr{Float64}, #unused#::Int32, grad_f::Ptr{Float64}, user_data::Ptr{Nothing})
    @ Ipopt ~/.julia/packages/Ipopt/M2QE8/src/C_wrapper.jl:49
 [39] IpoptSolve(prob::Ipopt.IpoptProblem)
    @ Ipopt ~/.julia/packages/Ipopt/M2QE8/src/C_wrapper.jl:433
 [40] optimize!(workspace::NonconvexIpopt.IpoptWorkspace{NonconvexCore.VecModel{NonconvexCore.Objective{NonconvexCore.var"#133#140"{Model{Vector{Any}}, NonconvexCore.Unflatten{Vector{Float64}, typeof(identity)}}, Base.RefValue{Float64}, Set{Symbol}}, NonconvexCore.VectorOfFunctions{Vector{NonconvexCore.EqConstraint}}, NonconvexCore.VectorOfFunctions{Vector{NonconvexCore.IneqConstraint{NonconvexCore.var"#137#144"{NonconvexCore.IneqConstraint{FunctionWrapper{var"#3#4", Int64}, Float64, Int64, Set{Symbol}}, NonconvexCore.Unflatten{Vector{Float64}, typeof(identity)}}, Float64, Int64, Set{Symbol}}}}, NonconvexCore.VectorOfFunctions{Vector{NonconvexCore.SDConstraint}}, Vector{Float64}, Vector{Float64}, Vector{Float64}}, Ipopt.IpoptProblem, Vector{Float64}, IpoptOptions{NamedTuple{(:tol, :hessian_approximation, :jac_c_constant, :jac_d_constant), Tuple{Float64, String, String, String}}}, Base.RefValue{Int64}})
    @ NonconvexIpopt ~/.julia/packages/NonconvexIpopt/jkxtN/src/ipopt.jl:52
 [41] #optimize#131
    @ ~/.julia/packages/NonconvexCore/uZAVo/src/models/vec_model.jl:63 [inlined]
 [42] optimize(::Model{Vector{Any}}, ::IpoptAlg, ::Vector{Float64}; kwargs::Base.Pairs{Symbol, IpoptOptions{NamedTuple{(:tol, :hessian_approximation, :jac_c_constant, :jac_d_constant), Tuple{Float64, String, String, String}}}, Tuple{Symbol}, NamedTuple{(:options,), Tuple{IpoptOptions{NamedTuple{(:tol, :hessian_approximation, :jac_c_constant, :jac_d_constant), Tuple{Float64, String, String, String}}}}}})
    @ NonconvexCore ~/.julia/packages/NonconvexCore/uZAVo/src/common.jl:233
 [43] top-level scope
    @ ~/Nextcloud/Dauphine/OptiForm/Julia/Nonconvex/testTrapz.jl:27
in expression starting at /home/mc/Nextcloud/Dauphine/OptiForm/Julia/Nonconvex/testTrapz.jl:27

```

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [June 14, 2022, 9:30am UTC](https://discourse.julialang.org/t/automatic-differentiation-and-trapz/82135/5 "2022-06-14T09:30:34Z")

</div>

The error is saying (unhelpfully) that `Trapz` isn’t compatible with automatic differentiation.

---

<div class="post-metadata">

**Author:** ![chupin](https://avatars.discourse-cdn.com/v4/letter/c/f07891/32.png) [@chupin](https://discourse.julialang.org/u/chupin)\
**Post date:** [June 14, 2022, 9:31am UTC](https://discourse.julialang.org/t/automatic-differentiation-and-trapz/82135/6 "2022-06-14T09:31:44Z")

</div>

Indeed 🙂 that’s why I posted this message.

---

<div class="post-metadata">

**Author:** ![mohamed82008](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mohamed82008/32/18171_2.png) [@mohamed82008](https://discourse.julialang.org/u/mohamed82008)\
**Post date:** [June 14, 2022, 12:22pm UTC](https://discourse.julialang.org/t/automatic-differentiation-and-trapz/82135/7 "2022-06-14T12:22:24Z")

</div>

How many decision variables are there? If not many, you can try to use ForwardDiff for the function. In the next release (in a day or 2), you can do `ff = forwarddiffy(f, x...)` to tell Nonconvex to use ForwardDiff for function `ff(x...)`.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [June 14, 2022, 1:04pm UTC](https://discourse.julialang.org/t/automatic-differentiation-and-trapz/82135/8 "2022-06-14T13:04:11Z")

</div>

> [@odow](#):
>
> The error is saying (unhelpfully) that `Trapz` isn’t compatible with automatic differentiation.

You could pretty easily define a manual differentiation rule (`rrule`) for `trapz` via ChainRules.jl. (This would make a good PR for the Trapz package too.)

Or you could use [Integrals.jl](https://github.com/SciML/Integrals.jl) instead of Trapz, since that provides differentiable integration routines.

---

<div class="post-metadata">

**Author:** ![antoine-levitt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/antoine-levitt/32/4008_2.png) [@antoine-levitt](https://discourse.julialang.org/u/antoine-levitt)\
**Post date:** [June 14, 2022, 2:43pm UTC](https://discourse.julialang.org/t/automatic-differentiation-and-trapz/82135/9 "2022-06-14T14:43:21Z")

</div>

Hi Maxime 🙂 To complete the responses here, the issue is that Zygote does not support mutation of arrays (which is anyway a pretty tricky thing to do in an AD engine). You can just precompute weights and replace `trapz(x)` with `sum(weights .* x)`: you’ll need to compute the weights if you want to write a `rrule` as @stevengj suggested anyway. Or just don’t use Trapz and use a simple uniform integration scheme.

---

<div class="post-metadata">

**Author:** ![chupin](https://avatars.discourse-cdn.com/v4/letter/c/f07891/32.png) [@chupin](https://discourse.julialang.org/u/chupin)\
**Post date:** [June 17, 2022, 9:46am UTC](https://discourse.julialang.org/t/automatic-differentiation-and-trapz/82135/10 "2022-06-17T09:46:44Z")

</div>

Thanks everyone. I’m going to test the different solutions.

---

<div class="post-metadata">

**Author:** ![mohamed82008](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mohamed82008/32/18171_2.png) [@mohamed82008](https://discourse.julialang.org/u/mohamed82008)\
**Post date:** [June 18, 2022, 5:53am UTC](https://discourse.julialang.org/t/automatic-differentiation-and-trapz/82135/11 "2022-06-18T05:53:18Z")

</div>

In Nonconvex 2.0, you can now use `abstractdiffy` to change the AD package used for a function or the entire model. See [Using other AD packages · Nonconvex.jl](https://julianonconvex.github.io/Nonconvex.jl/stable/gradients/other_ad/) for docs.