I am solving k-means using JuMP.

```
X = [
2.0 2.0;
2.5 2.1;
7.0 7.0;
2.2 2.3;
6.8 7.0;
7.2 7.5;]
k=2
m, n = size(X)
W = exp.(-Distances.pairwise(Distances.Euclidean(), X'))
I = Matrix(1.0 * LinearAlgebra.I, m, m)
model = JuMP.Model(SCS.Optimizer)
JuMP.unset_silent(model)
@JuMP.variable(model, Z[1:m, 1:m], PSD)
@JuMP.constraint(model, Z .>= 0)
@JuMP.constraint(model, sum(Z; dims=2) .== ones(m))
@JuMP.constraint(model, LinearAlgebra.tr(Z) == k)
@JuMP.objective(model, Min, LinearAlgebra.tr(W * (I - Z)))
JuMP.optimize!(model)
Z_value = JuMP.value.(Z)
```

Everything works great.

However, I want to impose constraints based on `Z`

value.

I wanted to add something like this:

```
JuMP.register(model, :proc_matrix, m * m, proc_matrix; autodiff = true)
@JuMP.variable(model, distances[1:m, 1:m] <= 0.1)
@JuMP.constraint(model, distances .== proc_matrix(Z...))
function proc_matrix(a...)
m = round(Int, floor(sqrt(length(a))))
M = reshape([a...], m, m)
Distances.pairwise(Distances.Euclidean(), M)
end
```

However, even this does not work.

Any suggestions?!

I also want to make `proc_matrix`

and the constraint more complicated but first I want to figure out how to process a variable matrix and return external output back to JuMP. Thank you!