# Zygote error - Mutating arrays is not supported

**URL:** <https://discourse.julialang.org/t/zygote-error-mutating-arrays-is-not-supported/89147>\
**Category:** New to Julia\
**Tags:** question, package, zygote\
**Created:** [October 23, 2022, 3:33pm UTC](https://discourse.julialang.org/t/zygote-error-mutating-arrays-is-not-supported/89147 "2022-10-23T15:33:59Z")\
**Posts on this page:** 13\
**Page:** 1

<div class="post-metadata">

**Author:** ![richinex](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/richinex/32/43805_2.png) [@richinex](https://discourse.julialang.org/u/richinex)\
**Post date:** [October 23, 2022, 3:33pm UTC](https://discourse.julialang.org/t/zygote-error-mutating-arrays-is-not-supported/89147/1 "2022-10-23T15:33:59Z")

</div>

Hi all, I am converting a function originally witten in python to Julia. but I have some problems working with the autodiff function provided by Zygote when I try to take the gradient of the function ( (I am able to do this seamlessly with Jax and pytorch)  
I understand that the lines where I fill in arrays P\_log and P\_norm are the problem but I am at a loss at what to do because P[kvals[i]:kvals[i + 1]-1] has a dynamic shape and sometimes produces a jagged array (i.e for each value of i, a different sized array might be produced - hence vcat is not an option).  
I will appreciate some some pointers on how to tackle this problem (Still trying to understand the inner workings of Julia)

```julia
function compute_total_obj(P::AbstractVector, 
    F::AbstractArray, 
    Z::AbstractArray, 
    Zerr_Re::AbstractArray,
    Zerr_Im::AbstractArray, 
    LB::AbstractVector, 
    UB::AbstractVector, 
    smf::AbstractVector, 
    func::Function, 
    num_params::Integer, 
    num_eis::Integer, 
    kvals::AbstractVector, 
    d2m::AbstractMatrix)

    # The lines causing the error
    P_log = zeros(num_params, num_eis)
    P_norm = zeros(num_params, num_eis)
    for i = 1:num_params
        P_log[i, :] .= (P[kvals[i]:kvals[i + 1]-1])

        P_norm[i, :] .= ((LB[kvals[i]:kvals[i + 1]-1] .+ (10 .^ P[kvals[i]:kvals[i + 1]-1])) ./ (1 .+ (10 .^ P[kvals[i]:kvals[i + 1]-1]) ./ UB[kvals[i]:kvals[i + 1]-1]))
    end
    smf_1 = ifelse.(isinf.(smf), 0.0, smf)
    chi_smf = sum(sum(((d2m * transpose(P_log)) .* (d2m * transpose(P_log))), dims = 1) .* smf_1)
    wrss_tot = compute_wrss.(eachcol(P_norm), eachcol(F), eachcol(Z), eachcol(Zerr_Re), eachcol(Zerr_Im), func)
    return (sum(wrss_tot) + chi_smf)
end

# Taking the gradient
    function g!(G, p)
        G .= Zygote.gradient(p -> compute_total_obj(p, F, Z, Zerr_Re, Zerr_Im, lb_vec, ub_vec, smf, func, num_params, num_eis, kvals, d2m), p)[1]
    end

# How kvals is formed (if any element of smf is Inf, a jagged array is produced
function get_kvals(smf::AbstractVector, num_eis::Integer)
    kvals = cumsum(insert!(ifelse.(isinf.(smf), 1, num_eis), 1, 1),)
    return kvals
end

# Error
# Output exceeds the size limit. Open the full output data in a text editor
# Mutating arrays is not supported -- called copyto!(SubArray{Float64, 1, Matrix{Float64}, Tuple{Int64, Base.Slice{Base.OneTo{Int64}}}, true}, ...)
# This error occurs when you ask Zygote to differentiate operations that change
# the elements of arrays in place (e.g. setting values with x .= ...)

# Possible fixes:
# - avoid mutating operations (preferred)
# - or read the documentation and solutions for this error
# https://fluxml.ai/Zygote.jl/latest/limitations

# Stacktrace:
# [1] error(s::String)
# @ Base .\error.jl:35
# [2] _throw_mutation_error(f::Function, args::SubArray{Float64, 1, Matrix{Float64}, Tuple{Int64, Base.Slice{Base.OneTo{Int64}}}, true})
# @ Zygote C:\Users\richinex\.julia\packages\Zygote\dABKa\src\lib\array.jl:68
# [3] (::Zygote.var"#389#390"{SubArray{Float64, 1, Matrix{Float64}, Tuple{Int64, Base.Slice{Base.OneTo{Int64}}}, true}})(#unused#::Nothing)
# @ Zygote C:\Users\richinex\.julia\packages\Zygote\dABKa\src\lib\array.jl:83
# [4] (::Zygote.var"#2474#back#391"{Zygote.var"#389#390"{SubArray{Float64, 1, Matrix{Float64}, Tuple{Int64, Base.Slice{Base.OneTo{Int64}}}, true}}})(Δ::Nothing)
# @ Zygote C:\Users\richinex\.julia\packages\ZygoteRules\AIbCs\src\adjoint.jl:67
# [5] Pullback
# @ .\broadcast.jl:871 [inlined]
# [6] Pullback
# @ .\broadcast.jl:868 [inlined]
# [7] Pullback
# @ .\broadcast.jl:864 [inlined]
# ...
# @ Optim C:\Users\richinex\.julia\packages\Optim\Zq1jM\src\multivariate\optimize\interface.jl:150
# [22] fit_deterministic(P::Vector{Float64}, freq::Vector{Float64}, Z::Matrix{ComplexF64}, bounds::Vector{Vector{Float64}}, smf::Vector{Float64}, func::Function; weight::Type)
# @ Main d:\julia_projects\error_calc_comparison\multi_2.ipynb:43
# [23] top-level scope
# @ d:\julia_projects\error_calc_comparison\multi_2.ipynb:1

```

---

<div class="post-metadata">

**Author:** ![ToucheSir](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/touchesir/32/14411_2.png) [@ToucheSir](https://discourse.julialang.org/u/ToucheSir)\
**Post date:** [October 23, 2022, 5:43pm UTC](https://discourse.julialang.org/t/zygote-error-mutating-arrays-is-not-supported/89147/2 "2022-10-23T17:43:28Z")

</div>

> [@richinex](#):
>
> ```julia
> P_log = zeros(num_params, num_eis)
> P_norm = zeros(num_params, num_eis)
> for i = 1:num_params
> P_log[i, :] .= (P[kvals[i]:kvals[i + 1]-1])
> 
> P_norm[i, :] .= ((LB[kvals[i]:kvals[i + 1]-1] .+ (10 .^ P[kvals[i]:kvals[i + 1]-1])) ./ (1 .+ (10 .^ P[kvals[i]:kvals[i + 1]-1]) ./ UB[kvals[i]:kvals[i + 1]-1]))
> end
> 
> ```

Given this wouldn’t work in JAX either (it might not error, but I’m pretty sure it’s at best undefined behaviour), this Julia code is presumably not a direct translation? If you can post a working Python example, we may be able to give some advice on a better one.

Edit: does the forward pass even work in Julia? This statement:

> [@richinex](#):
>
> P[kvals[i]:kvals[i + 1]-1] has a dynamic shape and sometimes produces a jagged array (i.e for each value of i, a different sized array might be produced - hence vcat is not an option

Is also true for the loop, because you can’t broadcast a length `M` array into a length `N > M` slice. So please do share a working Python example 🙂

---

<div class="post-metadata">

**Author:** ![richinex](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/richinex/32/43805_2.png) [@richinex](https://discourse.julialang.org/u/richinex)\
**Post date:** [October 23, 2022, 6:36pm UTC](https://discourse.julialang.org/t/zygote-error-mutating-arrays-is-not-supported/89147/3 "2022-10-23T18:36:51Z")

</div>

Hi ToucheSir, The code is pretty long and that is why I did not post it earlier. Here is the [link](https://colab.research.google.com/drive/1YZLfyYAIjUdOcR6TKQLUb4OAEmxMQutv?usp=sharing) to the complete code written in JAX which also contains the test data. I have just run it in google colab so I guess its working. Kindly go through and let me know how I can reproduce the compute\_total\_obj best in julia. Thanks.

---

<div class="post-metadata">

**Author:** ![ToucheSir](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/touchesir/32/14411_2.png) [@ToucheSir](https://discourse.julialang.org/u/ToucheSir)\
**Post date:** [October 23, 2022, 6:38pm UTC](https://discourse.julialang.org/t/zygote-error-mutating-arrays-is-not-supported/89147/4 "2022-10-23T18:38:29Z")

</div>

Doesn’t look like the notebook is publicly accessible, can you make it so, toss it in a GitHub gist or attach it here?

---

<div class="post-metadata">

**Author:** ![richinex](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/richinex/32/43805_2.png) [@richinex](https://discourse.julialang.org/u/richinex)\
**Post date:** [October 23, 2022, 7:00pm UTC](https://discourse.julialang.org/t/zygote-error-mutating-arrays-is-not-supported/89147/5 "2022-10-23T19:00:26Z")

</div>

> **[Google Colaboratory](https://colab.research.google.com/drive/1YZLfyYAIjUdOcR6TKQLUb4OAEmxMQutv?usp=sharing)**

---

<div class="post-metadata">

**Author:** ![richinex](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/richinex/32/43805_2.png) [@richinex](https://discourse.julialang.org/u/richinex)\
**Post date:** [October 23, 2022, 7:02pm UTC](https://discourse.julialang.org/t/zygote-error-mutating-arrays-is-not-supported/89147/6 "2022-10-23T19:02:35Z")

</div>

I’ve made it public.

Gist: [jax\_optim.ipynb · GitHub](https://gist.github.com/richinex/6a2ab066d040c0323aad1fe9dbb00f8e)

---

<div class="post-metadata">

**Author:** ![richinex](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/richinex/32/43805_2.png) [@richinex](https://discourse.julialang.org/u/richinex)\
**Post date:** [October 23, 2022, 8:43pm UTC](https://discourse.julialang.org/t/zygote-error-mutating-arrays-is-not-supported/89147/7 "2022-10-23T20:43:34Z")

</div>

To make comparison easier I have also attached my julia example notebook as a gist: [julia\_optim.ipynb · GitHub](https://gist.github.com/richinex/15b994727ebb0e1f8275576b4e17f0f9)

---

<div class="post-metadata">

**Author:** ![ToucheSir](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/touchesir/32/14411_2.png) [@ToucheSir](https://discourse.julialang.org/u/ToucheSir)\
**Post date:** [October 25, 2022, 5:17am UTC](https://discourse.julialang.org/t/zygote-error-mutating-arrays-is-not-supported/89147/8 "2022-10-25T05:17:13Z")

</div>

Thanks! Looking at a line like

```python
P_log = P_log.at[i, :].set(P[self.kvals[i]:self.kvals[i + 1]])

```

Is it correct to say that `len(P[self.kvals[i]:self.kvals[i + 1]]) == self.kvals[i + 1] - self.kvals[i] == self.num_eis`? If so, that should make the Julia code far easier to write. If not, you’ll have to explain to me how this code interacts with [🔪 JAX - The Sharp Bits 🔪 — JAX documentation](https://jax.readthedocs.io/en/latest/notebooks/Common_Gotchas_in_JAX.html?highlight=un-roll#out-of-bounds-indexing) 😄

---

<div class="post-metadata">

**Author:** ![richinex](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/richinex/32/43805_2.png) [@richinex](https://discourse.julialang.org/u/richinex)\
**Post date:** [October 25, 2022, 5:38am UTC](https://discourse.julialang.org/t/zygote-error-mutating-arrays-is-not-supported/89147/9 "2022-10-25T05:38:56Z")

</div>

Is it correct to say that `len(P[self.kvals[i]:self.kvals[i + 1]]) == self.kvals[i + 1] - self.kvals[i] == self.num_eis` ?  
This is only true when no value in `smf ` is Inf. when a value in smf is set to Inf, the intervals produced by the values in kval are not equal so I assign a scalar to P\_log.at[i,:] for instance something like P\_log.at[i,:].set(0.999). I found a workarround using julia’s ternary but it might be suboptimal. Maybe you get more hints on what I am trying to achieve and help to make it better. What I do is this: when the interval is just 1 value, I repeat that value to the size of num\_eis.

```julia
function compute_total_obj(P::AbstractVector, 
    F::AbstractArray, 
    Z::AbstractArray, 
    Zerr_Re::AbstractArray,
    Zerr_Im::AbstractArray, 
    LB::AbstractVector, 
    UB::AbstractVector, 
    smf::AbstractVector, 
    func::Function, 
    num_params::Integer, 
    num_eis::Integer, 
    kvals::AbstractVector, 
    d2m::AbstractMatrix)

    P_log = vcat(transpose.([length(P[kvals[i]:kvals[i + 1]-1]) == 1 ? repeat(P[kvals[i]:kvals[i + 1]-1], num_eis) : P[kvals[i]:kvals[i + 1]-1] for i = 1:num_params])...)
    P_norm = vcat(transpose.([length(P[kvals[i]:kvals[i + 1]-1]) == 1 ? repeat(((LB[kvals[i]:kvals[i + 1]-1] .+ (10 .^ P[kvals[i]:kvals[i + 1]-1])) ./ (1 .+ (10 .^ P[kvals[i]:kvals[i + 1]-1]) ./ UB[kvals[i]:kvals[i + 1]-1])), num_eis) : ((LB[kvals[i]:kvals[i + 1]-1] .+ (10 .^ P[kvals[i]:kvals[i + 1]-1])) ./ (1 .+ (10 .^ P[kvals[i]:kvals[i + 1]-1]) ./ UB[kvals[i]:kvals[i + 1]-1])) for i = 1:num_params])...)
    smf_1 = ifelse.(isinf.(smf), 0.0, smf)
    chi_smf = sum(sum(((d2m * transpose(P_log)) .* (d2m * transpose(P_log))), dims = 1) .* smf_1)
    wrss_tot = compute_wrss.(eachcol(P_norm), eachcol(F), eachcol(Z), eachcol(Zerr_Re), eachcol(Zerr_Im), func)
    return (sum(wrss_tot) + chi_smf)

```

---

<div class="post-metadata">

**Author:** ![richinex](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/richinex/32/43805_2.png) [@richinex](https://discourse.julialang.org/u/richinex)\
**Post date:** [October 25, 2022, 5:52am UTC](https://discourse.julialang.org/t/zygote-error-mutating-arrays-is-not-supported/89147/10 "2022-10-25T05:52:36Z")

</div>

The big idea is that I wanna be able to keep a parameter constant during optimization by setting it’s value in smf to Inf

---

<div class="post-metadata">

**Author:** ![ToucheSir](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/touchesir/32/14411_2.png) [@ToucheSir](https://discourse.julialang.org/u/ToucheSir)\
**Post date:** [October 25, 2022, 11:21pm UTC](https://discourse.julialang.org/t/zygote-error-mutating-arrays-is-not-supported/89147/11 "2022-10-25T23:21:52Z")

</div>

> [@richinex](#):
>
> when a value in smf is set to Inf, the intervals produced by the values in kval are not equal so I assign a scalar to P\_log.at[i,:] for instance something like P\_log.at[i,:].set(0.999).

I had to test this locally, but I think I understand now. JAX lets you do the equivalent of `array[:] = [x]`, so `len(P[self.kvals[i]:self.kvals[i + 1]])` may also be 1.

Your Julia implementation looks pretty reasonable. Here is my attempt at cleaning it up a bit:

```julia
P_log = reduce(hcat, [
    let istart=kvals[i], istop=kvals[i + 1] - 1, ps=@view P[istart:istop]
        istop - istart == 1 ? repeat(ps, num_eis) : ps
    end
    for i = 1:num_params
])
# Alternative way of doing this using map
# Not sure which is faster, feel free to try both
P_norm = reduce(hcat, map(1:num_params) do i
    istart, istop = kvals[i], kvals[i + 1] - 1
    # @view and @views reduces memory allocations when you're not differentiating this function with Zygote
    ps = @view P[istart:istop]
    ps_norm = @views (LB[istart:istop] .+ (10 .^ ps)) ./ (1 .+ (10 .^ ps) ./ UB[istart:istop])
    ub - lb == 1 ? repeat(ps, num_eis) : ps
end)

smf_1 = ifelse.(isinf.(smf), 0.0, smf)
chi_smf = sum(sum((d2m * transpose(P_log)).^2, dims = 1) .* smf_1)
wrss_tot = compute_wrss.(eachcol(P_norm), eachcol(F), eachcol(Z), eachcol(Zerr_Re), eachcol(Zerr_Im), func)
return sum(wrss_tot) + chi_smf

```

---

<div class="post-metadata">

**Author:** ![richinex](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/richinex/32/43805_2.png) [@richinex](https://discourse.julialang.org/u/richinex)\
**Post date:** [October 26, 2022, 7:14am UTC](https://discourse.julialang.org/t/zygote-error-mutating-arrays-is-not-supported/89147/12 "2022-10-26T07:14:01Z")

</div>

Thanks. Your clean up works but takes so much time (about 12 mins) to compile compared to using the one I posted earlier (using vcat) which takes less than 2 mins. I wonder why.

---

<div class="post-metadata">

**Author:** ![richinex](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/richinex/32/43805_2.png) [@richinex](https://discourse.julialang.org/u/richinex)\
**Post date:** [October 26, 2022, 7:29am UTC](https://discourse.julialang.org/t/zygote-error-mutating-arrays-is-not-supported/89147/13 "2022-10-26T07:29:10Z")

</div>

> [@ToucheSir](#):
>
> `ub - lb == 1 ? repeat(ps, num_eis) : ps`

It works now., I spotted the error. `ub - lb == 1 ? repeat(ps, num_eis) : ps` should be `istop - istart == 1 ? repeat(ps_norm, num_eis) : ps_norm`