Fix zygote constraint bug and update rosenbrock doc

docs/src/examples/rosenbrock.md

@@ -22,37 +22,45 @@ _p = [1.0, 100.0]
 f = OptimizationFunction(rosenbrock, Optimization.AutoForwardDiff())
 l1 = rosenbrock(x0, _p)
 prob = OptimizationProblem(f, x0, _p)
-
+```
 ## Optim.jl Solvers
 
-using OptimizationOptimJL
-
-# Start with some derivative-free optimizers
+### Start with some derivative-free optimizers
 
+```@example rosenbrock
+using OptimizationOptimJL
 sol = solve(prob, SimulatedAnnealing())
 prob = OptimizationProblem(f, x0, _p, lb = [-1.0, -1.0], ub = [0.8, 0.8])
 sol = solve(prob, SAMIN())
 
 l1 = rosenbrock(x0, _p)
 prob = OptimizationProblem(rosenbrock, x0, _p)
 sol = solve(prob, NelderMead())
+```
 
-# Now a gradient-based optimizer with forward-mode automatic differentiation
+### Now a gradient-based optimizer with forward-mode automatic differentiation
 
+```@example rosenbrock
 optf = OptimizationFunction(rosenbrock, Optimization.AutoForwardDiff())
 prob = OptimizationProblem(optf, x0, _p)
 sol = solve(prob, BFGS())
+```
 
-# Now a second order optimizer using Hessians generated by forward-mode automatic differentiation
+### Now a second order optimizer using Hessians generated by forward-mode automatic differentiation
 
+```@example rosenbrock
 sol = solve(prob, Newton())
+```
 
-# Now a second order Hessian-free optimizer
+### Now a second order Hessian-free optimizer
 
+```@example rosenbrock
 sol = solve(prob, Optim.KrylovTrustRegion())
+```
 
-# Now derivative-based optimizers with various constraints
+### Now derivative-based optimizers with various constraints
 
+```@example rosenbrock
 cons = (res, x, p) -> res .= [x[1]^2 + x[2]^2]
 optf = OptimizationFunction(rosenbrock, Optimization.AutoForwardDiff(); cons = cons)
 
@@ -68,24 +76,32 @@ sol = solve(prob, IPNewton())
 
 prob = OptimizationProblem(optf, x0, _p, lcons = [0.5], ucons = [0.5],
     lb = [-500.0, -500.0], ub = [50.0, 50.0])
-sol = solve(prob, IPNewton()) # Notice now that x[1]^2 + x[2]^2 ≈ 0.5:
-# cons(sol.u, _p) = 0.49999999999999994
+sol = solve(prob, IPNewton())
+
+# Notice now that x[1]^2 + x[2]^2 ≈ 0.5:
+cons(sol.u, _p)
+```
 
+```@example rosenbrock
 function con_c(res, x, p)
     res .= [x[1]^2 + x[2]^2]
 end
 
 optf = OptimizationFunction(rosenbrock, Optimization.AutoForwardDiff(); cons = con_c)
 prob = OptimizationProblem(optf, x0, _p, lcons = [-Inf], ucons = [0.25^2])
 sol = solve(prob, IPNewton()) # -Inf < cons_circ(sol.u, _p) = 0.25^2
+```
 
 ## Evolutionary.jl Solvers
 
+```@example rosenbrock
 using OptimizationEvolutionary
 sol = solve(prob, CMAES(μ = 40, λ = 100), abstol = 1e-15) # -Inf < cons_circ(sol.u, _p) = 0.25^2
+```
 
 ## IPOPT through OptimizationMOI
 
+```@example rosenbrock
 using OptimizationMOI, Ipopt
 
 function con2_c(res, x, p)
@@ -95,36 +111,46 @@ end
 optf = OptimizationFunction(rosenbrock, Optimization.AutoZygote(); cons = con2_c)
 prob = OptimizationProblem(optf, x0, _p, lcons = [-Inf, -Inf], ucons = [Inf, Inf])
 sol = solve(prob, Ipopt.Optimizer())
+```
 
-# Now let's switch over to OptimizationOptimisers with reverse-mode AD
+## Now let's switch over to OptimizationOptimisers with reverse-mode AD
 
+```@example rosenbrock
 using OptimizationOptimisers
 optf = OptimizationFunction(rosenbrock, Optimization.AutoZygote())
 prob = OptimizationProblem(optf, x0, _p)
 sol = solve(prob, Adam(0.05), maxiters = 1000, progress = false)
+```
 
 ## Try out CMAEvolutionStrategy.jl's evolutionary methods
 
+```@example rosenbrock
 using OptimizationCMAEvolutionStrategy
 sol = solve(prob, CMAEvolutionStrategyOpt())
+```
 
 ## Now try a few NLopt.jl solvers with symbolic differentiation via ModelingToolkit.jl
 
+```@example rosenbrock
 using OptimizationNLopt, ModelingToolkit
 optf = OptimizationFunction(rosenbrock, Optimization.AutoModelingToolkit())
 prob = OptimizationProblem(optf, x0, _p)
 
 sol = solve(prob, Opt(:LN_BOBYQA, 2))
 sol = solve(prob, Opt(:LD_LBFGS, 2))
+```
 
-## Add some box constraints and solve with a few NLopt.jl methods
+### Add some box constraints and solve with a few NLopt.jl methods
 
+```@example rosenbrock
 prob = OptimizationProblem(optf, x0, _p, lb = [-1.0, -1.0], ub = [0.8, 0.8])
 sol = solve(prob, Opt(:LD_LBFGS, 2))
 sol = solve(prob, Opt(:G_MLSL_LDS, 2), local_method = Opt(:LD_LBFGS, 2), maxiters = 10000) #a global optimizer with random starts of local optimization
+```
 
 ## BlackBoxOptim.jl Solvers
 
+```@example rosenbrock
 using OptimizationBBO
 prob = Optimization.OptimizationProblem(rosenbrock, x0, _p, lb = [-1.0, 0.2],
     ub = [0.8, 0.43])

ext/OptimizationZygoteExt.jl

@@ -123,13 +123,13 @@ function Optimization.instantiate_function(f, cache::Optimization.ReInitCache,
     if f.cons === nothing
         cons = nothing
     else
-        cons = (res, θ) -> f.cons(res, θ, cache.p)
-        cons_oop = (x) -> (_res = zeros(eltype(x), num_cons); cons(_res, x); _res)
+        cons =  (res, θ) -> f.cons(res, θ, cache.p)
+        cons_oop = (x) -> (_res = Zygote.Buffer(x, num_cons); cons(_res, x); copy(_res))
     end
 
     if cons !== nothing && f.cons_j === nothing
         cons_j = function (J, θ)
-            J .= Zygote.jacobian(cons_oop, θ)
+            J .= first(Zygote.jacobian(cons_oop, θ))
         end
     else
         cons_j = (J, θ) -> f.cons_j(J, θ, cache.p)