A small educational example:

```
using Zygote
using ForwardDiff
```

Let’s define a simple gaussian (sampled at integer coordinates)

```
x = -20:1.0:20
f2(xf) = 1000.0 .* exp.(-1.25 .* (x .- xf).^2)
```

And a simple center of mass computation to estimate the location of the peak from the samples

```
function cm(xc, α)
y = f2(xc)
# The next two lines are the issue
y = y .- 15.0;
y[y.<0] .= 0.0;
y = y.^α #just a random example, could be some other function
i = argmax(y)
idx = i-2:1:i+2
return sum(y[idx] .* x[idx]) / sum(y[idx])
end
```

And let’s say that I want to minimize the error that I make as function of alpha

```
function cost_cm(α)
r = cm.(-0.5:0.01:0.5, α) .- (-0.5:0.01:0.5)
return maximum(r)
end
```

Forward diff is fine

```
d_cost_cm(x) = ForwardDiff.derivative(cost_cm, x)
d_cost_cm(1.0)
-0.02776240571095623
```

Zygote is not happy

```
d_cost_cm_Zygote(x) = gradient(cost_cm, x)[1]
d_cost_cm_Zygote(1.0)
Mutating arrays is not supported ...
```

The issue are the lines:

```
y = y .- 15.0;
y[y.<0] .= 0.0;
```

Is there a way to replace them with something equivalent that would make Zygote happy?

(PS do we expect this issue to be solved by Diffractor.jl?)

Somehow related, if you use filter to create a new array instead of y[y.<0] .= 0.0 can you get the list of indexes to be used in sum(y_new_filtered_array .* x[idx]) / sum(y_new_filtered_array) ?