# When Zygote.jl meets Quadgk.jl: Mutating arrays is not supported -- called setindex!(::Vector{QuadGK.Segment})

**URL:** <https://discourse.julialang.org/t/when-zygote-jl-meets-quadgk-jl-mutating-arrays-is-not-supported-called-setindex-vector-quadgk-segment/84042>\
**Category:** General Usage\
**Tags:** zygote, quadgk\
**Created:** [July 11, 2022, 9:52am UTC](https://discourse.julialang.org/t/when-zygote-jl-meets-quadgk-jl-mutating-arrays-is-not-supported-called-setindex-vector-quadgk-segment/84042 "2022-07-11T09:52:01Z")\
**Posts on this page:** 7\
**Page:** 1

<div class="post-metadata">

**Author:** ![F-YF](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/f-yf/32/17363_2.png) [@F-YF](https://discourse.julialang.org/u/F-YF)\
**Post date:** [July 11, 2022, 9:52am UTC](https://discourse.julialang.org/t/when-zygote-jl-meets-quadgk-jl-mutating-arrays-is-not-supported-called-setindex-vector-quadgk-segment/84042/1 "2022-07-11T09:52:01Z")

</div>

I have a problem calculating the derivative of a function that involves an integral operation. Examples are as follows:

```julia
using QuadGK
using Zygote

function test()
    f(t)=quadgk(x->sin(x)*exp(t*cos(x)) , 0, pi, atol=1e-16)[1]
    @show f(1)
    g=gradient(x -> f(x), 1.0)[1]
    @show g
end
test()

```

The results are as follows:

```julia
f(1) = -4.428138526555854
ERROR: Mutating arrays is not supported -- called setindex!(::Vector{QuadGK.Segment{Float64, Float64, Float64}}, _...)
Stacktrace:
  [1] error(s::String)
    @ Base ./error.jl:33
  [2] (::Zygote.var"#441#442"{Vector{QuadGK.Segment{Float64, Float64, Float64}}})(#unused#::Nothing)
    @ Zygote ~/.julia/packages/Zygote/Y6SC4/src/lib/array.jl:71
  [3] (::Zygote.var"#2337#back#443"{Zygote.var"#441#442"{Vector{QuadGK.Segment{Float64, Float64, Float64}}}})(Δ::Nothing)
    @ Zygote ~/.julia/packages/ZygoteRules/AIbCs/src/adjoint.jl:67
  [4] Pullback
    @ ~/.julia/packages/DataStructures/izXYB/src/heaps/arrays_as_heaps.jl:45 [inlined]
  [5] (::typeof(∂(percolate_up!)))(Δ::Nothing)
    @ Zygote ~/.julia/packages/Zygote/Y6SC4/src/compiler/interface2.jl:0
  [6] Pullback
    @ ~/.julia/packages/DataStructures/izXYB/src/heaps/arrays_as_heaps.jl:48 [inlined]
  [7] Pullback
    @ ~/.julia/packages/DataStructures/izXYB/src/heaps/arrays_as_heaps.jl:73 [inlined]
  [8] (::typeof(∂(heappush!)))(Δ::Nothing)
    @ Zygote ~/.julia/packages/Zygote/Y6SC4/src/compiler/interface2.jl:0
  [9] Pullback
    @ ~/.julia/packages/QuadGK/ENhXl/src/adapt.jl:63 [inlined]
 [10] (::typeof(∂(adapt)))(Δ::Tuple{Float64, Nothing})
    @ Zygote ~/.julia/packages/Zygote/Y6SC4/src/compiler/interface2.jl:0
 [11] Pullback
    @ ~/.julia/packages/QuadGK/ENhXl/src/adapt.jl:35 [inlined]
 [12] (::typeof(∂(do_quadgk)))(Δ::Tuple{Float64, Nothing})
    @ Zygote ~/.julia/packages/Zygote/Y6SC4/src/compiler/interface2.jl:0
 [13] Pullback
    @ ~/.julia/packages/QuadGK/ENhXl/src/adapt.jl:179 [inlined]
 [14] (::typeof(∂(λ)))(Δ::Tuple{Float64, Nothing})
    @ Zygote ~/.julia/packages/Zygote/Y6SC4/src/compiler/interface2.jl:0
 [15] Pullback
    @ ~/.julia/packages/QuadGK/ENhXl/src/adapt.jl:113 [inlined]
 [16] (::typeof(∂(handle_infinities)))(Δ::Tuple{Float64, Nothing})
    @ Zygote ~/.julia/packages/Zygote/Y6SC4/src/compiler/interface2.jl:0
 [17] Pullback
    @ ~/.julia/packages/QuadGK/ENhXl/src/adapt.jl:177 [inlined]
 [18] (::typeof(∂(#quadgk#27)))(Δ::Tuple{Float64, Nothing})
    @ Zygote ~/.julia/packages/Zygote/Y6SC4/src/compiler/interface2.jl:0
 [19] #212
    @ ~/.julia/packages/Zygote/Y6SC4/src/lib/lib.jl:203 [inlined]
 [20] (::Zygote.var"#1750#back#214"{Zygote.var"#212#213"{Tuple{NTuple{9, Nothing}, Tuple{}}, typeof(∂(#quadgk#27))}})(Δ::Tuple{Float64, Nothing})
    @ Zygote ~/.julia/packages/ZygoteRules/AIbCs/src/adjoint.jl:67
 [21] Pullback
    @ ~/.julia/packages/QuadGK/ENhXl/src/adapt.jl:177 [inlined]
 [22] (::typeof(∂(quadgk##kw)))(Δ::Tuple{Float64, Nothing})
    @ Zygote ~/.julia/packages/Zygote/Y6SC4/src/compiler/interface2.jl:0
 [23] #212
    @ ~/.julia/packages/Zygote/Y6SC4/src/lib/lib.jl:203 [inlined]
 [24] (::Zygote.var"#1750#back#214"{Zygote.var"#212#213"{Tuple{Tuple{Nothing, Nothing, Nothing}, Tuple{Nothing, Nothing}}, typeof(∂(quadgk##kw))}})(Δ::Tuple{Float64, Nothing})
    @ Zygote ~/.julia/packages/ZygoteRules/AIbCs/src/adjoint.jl:67
 [25] Pullback
    @ ~/.julia/packages/QuadGK/ENhXl/src/adapt.jl:173 [inlined]
 [26] (::typeof(∂(#quadgk#26)))(Δ::Tuple{Float64, Nothing})
    @ Zygote ~/.julia/packages/Zygote/Y6SC4/src/compiler/interface2.jl:0
 [27] #212
    @ ~/.julia/packages/Zygote/Y6SC4/src/lib/lib.jl:203 [inlined]
 [28] #1750#back
    @ ~/.julia/packages/ZygoteRules/AIbCs/src/adjoint.jl:67 [inlined]
 [29] Pullback
    @ ~/.julia/packages/QuadGK/ENhXl/src/adapt.jl:173 [inlined]
 [30] (::typeof(∂(quadgk##kw)))(Δ::Tuple{Float64, Nothing})
    @ Zygote ~/.julia/packages/Zygote/Y6SC4/src/compiler/interface2.jl:0
 [31] Pullback
    @ ./REPL[3]:2 [inlined]
 [32] (::typeof(∂(f)))(Δ::Float64)
    @ Zygote ~/.julia/packages/Zygote/Y6SC4/src/compiler/interface2.jl:0
 [33] Pullback
    @ ./REPL[3]:4 [inlined]
 [34] (::Zygote.var"#56#57"{typeof(∂(#2))})(Δ::Float64)
    @ Zygote ~/.julia/packages/Zygote/Y6SC4/src/compiler/interface.jl:41
 [35] gradient(f::Function, args::Float64)
    @ Zygote ~/.julia/packages/Zygote/Y6SC4/src/compiler/interface.jl:76
 [36] test()
    @ Main ./REPL[3]:4
 [37] top-level scope
    @ REPL[4]:1

```

Zygote.buffer() doesn’t seem to solve this problem. I want to consult what good idea does everybody have to be able to solve the problem.

---

<div class="post-metadata">

**Author:** ![bertschi](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bertschi/32/33462_2.png) [@bertschi](https://discourse.julialang.org/u/bertschi)\
**Post date:** [July 11, 2022, 6:24pm UTC](https://discourse.julialang.org/t/when-zygote-jl-meets-quadgk-jl-mutating-arrays-is-not-supported-called-setindex-vector-quadgk-segment/84042/2 "2022-07-11T18:24:32Z")

</div>

In your example `f(x)` is not actually a function of `x`, but rather a constant as your range of integration is fixed:

\int\_0^\pi e^x \cdot \sin x \cdot \cos x \; dx \approx -4.428

Thus, you don’t need to take any derivative.  
In case, your integration range depends on `x` I would suggest adding an adjoint for `quadgk` which makes use of the fundamental theorem of calculus. This should be way more efficient than trying to autodiff through a numeric integral.

---

<div class="post-metadata">

**Author:** ![tbeason](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tbeason/32/15898_2.png) [@tbeason](https://discourse.julialang.org/u/tbeason)\
**Post date:** [July 11, 2022, 6:47pm UTC](https://discourse.julialang.org/t/when-zygote-jl-meets-quadgk-jl-mutating-arrays-is-not-supported-called-setindex-vector-quadgk-segment/84042/3 "2022-07-11T18:47:53Z")

</div>

QuadGK is not differentiable.

---

<div class="post-metadata">

**Author:** ![jlperla](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jlperla/32/34332_2.png) [@jlperla](https://discourse.julialang.org/u/jlperla)\
**Post date:** [July 11, 2022, 7:39pm UTC](https://discourse.julialang.org/t/when-zygote-jl-meets-quadgk-jl-mutating-arrays-is-not-supported-called-setindex-vector-quadgk-segment/84042/4 "2022-07-11T19:39:19Z")

</div>

I think this is what [GitHub - SciML/Integrals.jl: A common interface for quadrature and numerical integration for the SciML scientific machine learning organization](https://github.com/SciML/Integrals.jl) etc. is intended for, but not sure the exact status of available algorithms. Looks like QuadGK is available, though? [Integral Solver Algorithms · Integrals.jl](https://integrals.sciml.ai/dev/solvers/IntegralSolvers/)

Regadless, you don’t want to differentiate the quadrature directly and need an adjoint rule implementing the calculus at a higher level.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [July 11, 2022, 11:45pm UTC](https://discourse.julialang.org/t/when-zygote-jl-meets-quadgk-jl-mutating-arrays-is-not-supported-called-setindex-vector-quadgk-segment/84042/5 "2022-07-11T23:45:30Z")

</div>

> [@jlperla](#):
>
> I think this is what [GitHub - SciML/Integrals.jl: A common interface for quadrature and numerical integration for the SciML scientific machine learning organization](https://github.com/SciML/Integrals.jl) etc. is intended for, but not sure the exact status of available algorithms. Looks like QuadGK is available, though? [Integral Solver Algorithms · Integrals.jl](https://integrals.sciml.ai/dev/solvers/IntegralSolvers/)

yeah it defines the differentiation overloads. See:

> <https://github.com/SciML/Integrals.jl/blob/master/test/derivative_tests.jl#L7-L29>

It doesn’t correct differentiation of the integral limits with forward mode right now, but it does do them with reverse. And the derivative w.r.t. parameters is handled in both cases.

---

<div class="post-metadata">

**Author:** ![F-YF](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/f-yf/32/17363_2.png) [@F-YF](https://discourse.julialang.org/u/F-YF)\
**Post date:** [July 12, 2022, 1:16am UTC](https://discourse.julialang.org/t/when-zygote-jl-meets-quadgk-jl-mutating-arrays-is-not-supported-called-setindex-vector-quadgk-segment/84042/6 "2022-07-12T01:16:36Z")

</div>

> [@bertschi](#):
>
> In your example `f(x)` is not actually a function of `x`, but rather a constant as your range of integration is fixed:
> 
> \int\_0^\pi e^x \cdot \sin x \cdot \cos x ; dx \approx -4.428
> 
> ∫π0ex⋅sinx⋅cosxdx≈−4.428 \int\_0^\pi e^x \cdot \sin x \cdot \cos x ; dx \approx -4.428
> 
> Thus, you don’t need to take any derivative.

I’m terribly sorry that it is my mistake. I think my function is going to look something like this,

```julia
using QuadGK
using Zygote

function test()
    f(t)=quadgk(x->sin(x)*exp(t*cos(x)) , 0, pi, atol=1e-16)[1]
    @show f(1)
    g=gradient(x -> f(x), 1.0)[1]
    @show g
end
test()

```

- What I really need is something like this

---

<div class="post-metadata">

**Author:** ![F-YF](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/f-yf/32/17363_2.png) [@F-YF](https://discourse.julialang.org/u/F-YF)\
**Post date:** [July 12, 2022, 5:14am UTC](https://discourse.julialang.org/t/when-zygote-jl-meets-quadgk-jl-mutating-arrays-is-not-supported-called-setindex-vector-quadgk-segment/84042/7 "2022-07-12T05:14:15Z")

</div>

> [@jlperla](#):
>
> I think this is what [GitHub - SciML/Integrals.jl: A common interface for quadrature and numerical integration for the SciML scientific machine learning organization](https://github.com/SciML/Integrals.jl) etc. is intended for, but not sure the exact status of available algorithms. Looks like QuadGK is available, though? [Integral Solver Algorithms · Integrals.jl](https://integrals.sciml.ai/dev/solvers/IntegralSolvers/)
> 
> Regadless, you don’t want to differentiate the quadrature directly and need an adjoint rule implementing the calculus at a higher level.

It’s great that [Integrals.jl](https://integrals.sciml.ai/dev/solvers/IntegralSolvers/) solved my problem. Thank you very much!