diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index 36f262b35..591f9917d 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -34,6 +34,7 @@ jobs: - OptimizationOptimisers - OptimizationPRIMA - OptimizationQuadDIRECT + - OptimizationSciPy - OptimizationSpeedMapping - OptimizationPolyalgorithms - OptimizationNLPModels diff --git a/lib/OptimizationSciPy/CondaPkg.toml b/lib/OptimizationSciPy/CondaPkg.toml new file mode 100644 index 000000000..7644a1f81 --- /dev/null +++ b/lib/OptimizationSciPy/CondaPkg.toml @@ -0,0 +1,3 @@ +[deps] +scipy = "" +numpy = "" \ No newline at end of file diff --git a/lib/OptimizationSciPy/LICENSE b/lib/OptimizationSciPy/LICENSE new file mode 100644 index 000000000..4647d51e1 --- /dev/null +++ b/lib/OptimizationSciPy/LICENSE @@ -0,0 +1,21 @@ +MIT License + +Copyright (c) 2023 Vaibhav Dixit and contributors + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. \ No newline at end of file diff --git a/lib/OptimizationSciPy/Project.toml b/lib/OptimizationSciPy/Project.toml new file mode 100644 index 000000000..ed9c72d71 --- /dev/null +++ b/lib/OptimizationSciPy/Project.toml @@ -0,0 +1,26 @@ +name = "OptimizationSciPy" +uuid = "cce07bd8-c79b-4b00-aee8-8db9cce22837" +authors = ["Aditya Pandey and contributors"] +version = "0.4.1" + +[deps] +Optimization = "7f7a1694-90dd-40f0-9382-eb1efda571ba" +PythonCall = "6099a3de-0909-46bc-b1f4-468b9a2dfc0d" +Reexport = "189a3867-3050-52da-a836-e630ba90ab69" + +[compat] +Optimization = "4" +PythonCall = "0.9" +Reexport = "1.2" +julia = "1" + +[extras] +ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" +ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78" +Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +ReverseDiff = "37e2e3b7-166d-5795-8a7a-e32c996b4267" +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" +Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f" + +[targets] +test = ["ForwardDiff", "ModelingToolkit", "Random", "ReverseDiff", "Test", "Zygote"] \ No newline at end of file diff --git a/lib/OptimizationSciPy/src/OptimizationSciPy.jl b/lib/OptimizationSciPy/src/OptimizationSciPy.jl new file mode 100644 index 000000000..e153a8bb9 --- /dev/null +++ b/lib/OptimizationSciPy/src/OptimizationSciPy.jl @@ -0,0 +1,1484 @@ +#This file lets you drive SciPy optimizers through SciML's Optimization.jl API. +module OptimizationSciPy + +using Reexport +@reexport using Optimization +using Optimization.SciMLBase +using PythonCall + +# We keep a handle to the actual Python SciPy module here. +const scipy = PythonCall.pynew() + +function __init__() + PythonCall.pycopy!(scipy, pyimport("scipy")) +end + +# Make sure whatever we got back is a plain Julia Vector{T}. +function ensure_julia_array(x, ::Type{T}=Float64) where T + x isa Vector{T} && return x + return convert(Vector{T}, x isa Py ? pyconvert(Vector, x) : x) +end + +# Pull a human-readable message out of the SciPy result object. +function safe_get_message(result) + pyhasattr(result, "message") || return "Optimization completed" + msg = result.message + if pyisinstance(msg, pybuiltins.str) + return pyconvert(String, msg) + end + if pyisinstance(msg, pybuiltins.list) || pyisinstance(msg, pybuiltins.tuple) + return join(pyconvert(Vector{String}, msg), ", ") + end + return string(pytypeof(msg)) +end + +# Squash any kind of numeric object down to a Julia Float64. +function safe_to_float(x) + x isa Float64 && return x + x isa Number && return Float64(x) + + if x isa Py + if pyhasattr(x, "item") + v = pyconvert(Float64, x.item(), nothing) + v !== nothing && return v + end + v = pyconvert(Float64, x, nothing) + v !== nothing && return v + end + + error("Cannot convert object to Float64: $(typeof(x))") +end + +# Gather timing / iteration counts and wrap them in OptimizationStats. +function extract_stats(result, time_elapsed) + stats_dict = Dict{Symbol, Any}( + :iterations => 0, + :time => time_elapsed, + :fevals => 0, + :gevals => 0, + :hevals => 0 + ) + if pyhasattr(result, "nit") && !pyis(result.nit, pybuiltins.None) + stats_dict[:iterations] = pyconvert(Int, result.nit) + end + if pyhasattr(result, "nfev") && !pyis(result.nfev, pybuiltins.None) + stats_dict[:fevals] = pyconvert(Int, result.nfev) + end + if pyhasattr(result, "njev") && !pyis(result.njev, pybuiltins.None) + stats_dict[:gevals] = pyconvert(Int, result.njev) + elseif pyhasattr(result, "ngrad") && !pyis(result.ngrad, pybuiltins.None) + stats_dict[:gevals] = pyconvert(Int, result.ngrad) + end + if pyhasattr(result, "nhev") && !pyis(result.nhev, pybuiltins.None) + stats_dict[:hevals] = pyconvert(Int, result.nhev) + end + return Optimization.OptimizationStats(; stats_dict...) +end + +# Map SciPy status integers onto SciML ReturnCode symbols. +function scipy_status_to_retcode(status::Int, success::Bool) + if success + return SciMLBase.ReturnCode.Success + end + return if status == 0 + SciMLBase.ReturnCode.Success + elseif status == 1 + SciMLBase.ReturnCode.MaxIters + elseif status == 2 + SciMLBase.ReturnCode.Infeasible + elseif status == 3 + SciMLBase.ReturnCode.Unstable + elseif status == 4 + SciMLBase.ReturnCode.Terminated + elseif status == 9 + SciMLBase.ReturnCode.MaxIters + else + SciMLBase.ReturnCode.Failure + end +end + +# Tiny structs that tag which SciPy algorithm the user picked. +abstract type ScipyOptimizer end + +struct ScipyMinimize <: ScipyOptimizer + method::String + function ScipyMinimize(method::String) + valid_methods = ["Nelder-Mead", "Powell", "CG", "BFGS", "Newton-CG", + "L-BFGS-B", "TNC", "COBYLA", "COBYQA", "SLSQP", + "trust-constr", "dogleg", "trust-ncg", "trust-krylov", + "trust-exact"] + if !(method in valid_methods) + throw(ArgumentError("Invalid method: $method. Valid methods are: $(join(valid_methods, ", "))")) + end + new(method) + end +end +ScipyMinimize() = ScipyMinimize("BFGS") + +ScipyNelderMead() = ScipyMinimize("Nelder-Mead") +ScipyPowell() = ScipyMinimize("Powell") +ScipyCG() = ScipyMinimize("CG") +ScipyBFGS() = ScipyMinimize("BFGS") +ScipyNewtonCG() = ScipyMinimize("Newton-CG") +ScipyLBFGSB() = ScipyMinimize("L-BFGS-B") +ScipyTNC() = ScipyMinimize("TNC") +ScipyCOBYLA() = ScipyMinimize("COBYLA") +ScipyCOBYQA() = ScipyMinimize("COBYQA") +ScipySLSQP() = ScipyMinimize("SLSQP") +ScipyTrustConstr() = ScipyMinimize("trust-constr") +ScipyDogleg() = ScipyMinimize("dogleg") +ScipyTrustNCG() = ScipyMinimize("trust-ncg") +ScipyTrustKrylov() = ScipyMinimize("trust-krylov") +ScipyTrustExact() = ScipyMinimize("trust-exact") + +struct ScipyMinimizeScalar <: ScipyOptimizer + method::String + function ScipyMinimizeScalar(method::String="brent") + valid_methods = ["brent", "bounded", "golden"] + if !(method in valid_methods) + throw(ArgumentError("Invalid method: $method. Valid methods are: $(join(valid_methods, ", "))")) + end + new(method) + end +end + +ScipyBrent() = ScipyMinimizeScalar("brent") +ScipyBounded() = ScipyMinimizeScalar("bounded") +ScipyGolden() = ScipyMinimizeScalar("golden") + +struct ScipyLeastSquares <: ScipyOptimizer + method::String + loss::String + function ScipyLeastSquares(; method::String="trf", loss::String="linear") + valid_methods = ["trf", "dogbox", "lm"] + valid_losses = ["linear", "soft_l1", "huber", "cauchy", "arctan"] + if !(method in valid_methods) + throw(ArgumentError("Invalid method: $method. Valid methods are: $(join(valid_methods, ", "))")) + end + if !(loss in valid_losses) + throw(ArgumentError("Invalid loss: $loss. Valid loss functions are: $(join(valid_losses, ", "))")) + end + new(method, loss) + end +end + +ScipyLeastSquaresTRF() = ScipyLeastSquares(method="trf") +ScipyLeastSquaresDogbox() = ScipyLeastSquares(method="dogbox") +ScipyLeastSquaresLM() = ScipyLeastSquares(method="lm") + +struct ScipyRootScalar <: ScipyOptimizer + method::String + function ScipyRootScalar(method::String="brentq") + valid_methods = ["brentq", "brenth", "bisect", "ridder", "newton", "secant", "halley", "toms748"] + if !(method in valid_methods) + throw(ArgumentError("Invalid method: $method. Valid methods are: $(join(valid_methods, ", "))")) + end + new(method) + end +end + +struct ScipyRoot <: ScipyOptimizer + method::String + function ScipyRoot(method::String="hybr") + valid_methods = ["hybr", "lm", "broyden1", "broyden2", "anderson", + "linearmixing", "diagbroyden", "excitingmixing", + "krylov", "df-sane"] + if !(method in valid_methods) + throw(ArgumentError("Invalid method: $method. Valid methods are: $(join(valid_methods, ", "))")) + end + new(method) + end +end + +struct ScipyLinprog <: ScipyOptimizer + method::String + function ScipyLinprog(method::String="highs") + valid_methods = ["highs", "highs-ds", "highs-ipm", "interior-point", + "revised simplex", "simplex"] + if !(method in valid_methods) + throw(ArgumentError("Invalid method: $method. Valid methods are: $(join(valid_methods, ", "))")) + end + new(method) + end +end + +struct ScipyMilp <: ScipyOptimizer end +struct ScipyDifferentialEvolution <: ScipyOptimizer end +struct ScipyBasinhopping <: ScipyOptimizer end +struct ScipyDualAnnealing <: ScipyOptimizer end +struct ScipyShgo <: ScipyOptimizer end +struct ScipyDirect <: ScipyOptimizer end +struct ScipyBrute <: ScipyOptimizer end + +for opt_type in [:ScipyMinimize, :ScipyDifferentialEvolution, :ScipyBasinhopping, + :ScipyDualAnnealing, :ScipyShgo, :ScipyDirect, :ScipyBrute, + :ScipyLinprog, :ScipyMilp] + @eval begin + SciMLBase.allowsbounds(::$opt_type) = true + SciMLBase.supports_opt_cache_interface(::$opt_type) = true + end +end + +for opt_type in [:ScipyMinimizeScalar, :ScipyRootScalar, :ScipyLeastSquares] + @eval begin + SciMLBase.supports_opt_cache_interface(::$opt_type) = true + end +end + +SciMLBase.supports_opt_cache_interface(::ScipyRoot) = true + +function SciMLBase.requiresgradient(opt::ScipyMinimize) + gradient_free = ["Nelder-Mead", "Powell", "COBYLA", "COBYQA"] + return !(opt.method in gradient_free) +end + +for opt_type in [:ScipyDifferentialEvolution, :ScipyBasinhopping, + :ScipyDualAnnealing, :ScipyShgo, :ScipyDirect, :ScipyBrute, + :ScipyMinimizeScalar, :ScipyLeastSquares, :ScipyRootScalar, + :ScipyRoot, :ScipyLinprog, :ScipyMilp] + @eval SciMLBase.requiresgradient(::$opt_type) = false +end + +function SciMLBase.requireshessian(opt::ScipyMinimize) + hessian_methods = ["Newton-CG", "dogleg", "trust-ncg", "trust-exact", "trust-krylov"] + return opt.method in hessian_methods +end + +function SciMLBase.requireshessian(opt::ScipyRootScalar) + return opt.method == "halley" +end + +function SciMLBase.allowsconstraints(opt::ScipyMinimize) + return opt.method in ["SLSQP", "trust-constr", "COBYLA", "COBYQA"] +end + +function SciMLBase.requiresconsjac(opt::ScipyMinimize) + return opt.method in ["SLSQP", "trust-constr"] +end + +SciMLBase.allowsconstraints(::ScipyShgo) = true +SciMLBase.allowsconstraints(::ScipyLinprog) = true +SciMLBase.allowsconstraints(::ScipyMilp) = true + +function SciMLBase.allowsbounds(opt::ScipyMinimizeScalar) + return opt.method == "bounded" +end + +function SciMLBase.allowsbounds(opt::ScipyLeastSquares) + return opt.method in ["trf", "dogbox"] +end + +function SciMLBase.allowsbounds(opt::ScipyRootScalar) + return opt.method in ["brentq", "brenth", "bisect", "ridder"] +end + +SciMLBase.allowsbounds(::ScipyRoot) = false + +function SciMLBase.__init(prob::SciMLBase.OptimizationProblem, opt::ScipyOptimizer; + cons_tol = 1e-6, + callback = (args...) -> (false), + progress = false, + kwargs...) + requires_bounds = opt isa Union{ScipyDifferentialEvolution, ScipyDirect, ScipyDualAnnealing, ScipyBrute} + if requires_bounds && (isnothing(prob.lb) || isnothing(prob.ub)) + throw(SciMLBase.IncompatibleOptimizerError("$(typeof(opt)) requires bounds")) + end + if opt isa ScipyMinimizeScalar && length(prob.u0) != 1 + throw(ArgumentError("ScipyMinimizeScalar requires exactly 1 variable, got $(length(prob.u0)). Use ScipyMinimize for multivariate problems.")) + end + if opt isa ScipyRootScalar && length(prob.u0) != 1 + throw(ArgumentError("ScipyRootScalar requires exactly 1 variable, got $(length(prob.u0)). Use ScipyRoot for multivariate problems.")) + end + if opt isa ScipyMinimizeScalar && opt.method == "bounded" + if isnothing(prob.lb) || isnothing(prob.ub) + throw(ArgumentError("ScipyMinimizeScalar with method='bounded' requires bounds")) + end + end + if opt isa ScipyRootScalar && opt.method in ["brentq", "brenth", "bisect", "ridder"] + if isnothing(prob.lb) || isnothing(prob.ub) + throw(ArgumentError("ScipyRootScalar with method='$(opt.method)' requires bracket (bounds)")) + end + end + if !isnothing(prob.lb) && !isnothing(prob.ub) + @assert length(prob.lb) == length(prob.ub) "Bounds must have the same length" + @assert all(prob.lb .<= prob.ub) "Lower bounds must be less than or equal to upper bounds" + end + return OptimizationCache(prob, opt; cons_tol, callback, progress, kwargs...) +end + +function SciMLBase.__solve(cache::OptimizationCache{F,RC,LB,UB,LC,UC,S,O,D,P,C}) where + {F,RC,LB,UB,LC,UC,S,O<:ScipyMinimize,D,P,C} + local cons_cache = nothing + if !isnothing(cache.f.cons) && !isnothing(cache.lcons) + cons_cache = zeros(eltype(cache.u0), length(cache.lcons)) + end + _loss = _create_loss(cache) + maxiters = Optimization._check_and_convert_maxiters(cache.solver_args.maxiters) + abstol = cache.solver_args.abstol + reltol = cache.solver_args.reltol + options = Dict{String, Any}() + if cache.opt.method == "trust-constr" + options["initial_tr_radius"] = 1.0 + options["verbose"] = 0 + options["finite_diff_rel_step"] = 1e-8 + options["gtol"] = 1e-10 + options["maxiter"] = 50000 + elseif cache.opt.method in ["dogleg", "trust-ncg", "trust-krylov", "trust-exact"] + options["gtol"] = 1e-10 + options["maxiter"] = 50000 + end + if !isnothing(maxiters) + options["maxiter"] = maxiters + end + if !isnothing(abstol) + if cache.opt.method in ["Nelder-Mead", "Powell"] + options["xatol"] = abstol + elseif cache.opt.method in ["L-BFGS-B", "TNC", "SLSQP", "trust-constr"] + options["ftol"] = abstol + elseif cache.opt.method == "COBYQA" + options["feasibility_tol"] = abstol + end + end + if !isnothing(reltol) + if cache.opt.method in ["CG", "BFGS", "Newton-CG", "L-BFGS-B", "TNC", "SLSQP", "trust-constr"] + options["gtol"] = reltol + end + end + _merge_solver_kwargs!(options, cache.solver_args) + jac = nothing + if SciMLBase.requiresgradient(cache.opt) && !isnothing(cache.f.grad) + _grad = function(θ) + θ_julia = ensure_julia_array(θ, eltype(cache.u0)) + grad = zeros(eltype(cache.u0), length(θ_julia)) + cache.f.grad(grad, θ_julia, cache.p) + return cache.sense === Optimization.MaxSense ? -grad : grad + end + jac = _grad + end + hess = nothing + if SciMLBase.requireshessian(cache.opt) + if !isnothing(cache.f.hess) + _hess = function(θ) + θ_julia = ensure_julia_array(θ, eltype(cache.u0)) + H = zeros(eltype(cache.u0), length(θ_julia), length(θ_julia)) + cache.f.hess(H, θ_julia, cache.p) + return cache.sense === Optimization.MaxSense ? -H : H + end + hess = _hess + else + if cache.opt.method in ["trust-constr", "dogleg", "trust-ncg", "trust-krylov", "trust-exact"] + options["hess"] = "BFGS" + else + throw(ArgumentError("Method $(cache.opt.method) requires Hessian but none was provided")) + end + end + end + bounds = nothing + if !isnothing(cache.lb) && !isnothing(cache.ub) + if cache.opt.method in ["L-BFGS-B", "TNC", "SLSQP", "trust-constr", "COBYLA", "COBYQA"] + bounds = scipy.optimize.Bounds(cache.lb, cache.ub) + end + end + constraints = pylist([]) + if SciMLBase.allowsconstraints(cache.opt) + if !isnothing(cache.f.cons) && !isnothing(cons_cache) + lcons = cache.lcons + ucons = cache.ucons + _cons_func = function(θ) + θ_julia = ensure_julia_array(θ, eltype(cache.u0)) + cons_cache .= zero(eltype(cons_cache)) + if hasmethod(cache.f.cons, Tuple{typeof(cons_cache), typeof(θ_julia), typeof(cache.p)}) + cache.f.cons(cons_cache, θ_julia, cache.p) + else + cache.f.cons(cons_cache, θ_julia) + end + return cons_cache + end + cons_jac = "2-point" + if SciMLBase.requiresconsjac(cache.opt) && !isnothing(cache.f.cons_j) + cons_j_cache = zeros(eltype(cache.u0), length(lcons), length(cache.u0)) + _cons_jac = function(θ) + θ_julia = ensure_julia_array(θ, eltype(cache.u0)) + if hasmethod(cache.f.cons_j, Tuple{typeof(cons_j_cache), typeof(θ_julia), typeof(cache.p)}) + cache.f.cons_j(cons_j_cache, θ_julia, cache.p) + else + cache.f.cons_j(cons_j_cache, θ_julia) + end + return cons_j_cache + end + cons_jac = _cons_jac + end + # user-controlled NonlinearConstraint extras + keep_feasible_flag = get(cache.solver_args, :keep_feasible, false) + jac_sparsity = get(cache.solver_args, :jac_sparsity, nothing) + nlc = scipy.optimize.NonlinearConstraint( + _cons_func, + lcons, + ucons; + jac = cons_jac, + keep_feasible = keep_feasible_flag, + finite_diff_rel_step = get(cache.solver_args, :cons_tol, 1e-8), + finite_diff_jac_sparsity = jac_sparsity, + ) + constraints = pylist([nlc]) + end + elseif !isnothing(cache.f.cons) + throw(ArgumentError("Method $(cache.opt.method) does not support constraints. Use SLSQP, trust-constr, COBYLA, or COBYQA instead.")) + end + # allow users to specify a Hessian update strategy (e.g. "BFGS", "SR1") + if cache.opt.method == "trust-constr" + hess_update = get(cache.solver_args, :hess_update, nothing) + if hess_update !== nothing + hess = hess_update + end + end + t0 = time() + result = nothing + try + result = scipy.optimize.minimize( + _loss, + cache.u0, + method = cache.opt.method, + jac = jac, + hess = hess, + bounds = bounds, + constraints = constraints, + options = pydict(options) + ) + catch e + if e isa PythonCall.Core.PyException + py_msg = sprint(showerror, e) + if occursin("Optimization halted by callback", py_msg) + throw(ErrorException("Optimization halted by callback")) + elseif occursin("Optimization halted: time limit exceeded", py_msg) + throw(ErrorException("Optimization halted: time limit exceeded")) + else + throw(ErrorException("SciPy optimization failed: $py_msg")) + end + else + rethrow(e) + end + end + if isnothing(result) + throw(ErrorException("Optimization failed to return a result")) + end + t1 = time() + if pyis(result.x, pybuiltins.None) + minimizer = fill(NaN, length(cache.u0)) + else + minimizer = pyconvert(Vector{eltype(cache.u0)}, result.x) + end + minimum = pyis(result.fun, pybuiltins.None) ? NaN : safe_to_float(result.fun) + py_success = pyconvert(Bool, pybool(result.success)) + py_message = safe_get_message(result) + status = 0 + if pyhasattr(result, "status") + try + status = pyconvert(Int, result.status) + catch + end + end + if cache.sense === Optimization.MaxSense + minimum = -minimum + end + retcode = scipy_status_to_retcode(status, py_success) + if retcode != SciMLBase.ReturnCode.Success + @debug "ScipyMinimize convergence: $(py_message)" + end + stats = extract_stats(result, t1 - t0) + return SciMLBase.build_solution(cache, cache.opt, minimizer, minimum; + original = result, + retcode = retcode, + stats = stats) +end + +function SciMLBase.__solve(cache::OptimizationCache{F,RC,LB,UB,LC,UC,S,O,D,P,C}) where + {F,RC,LB,UB,LC,UC,S,O<:ScipyMinimizeScalar,D,P,C} + maxtime = get(cache.solver_args, :maxtime, nothing) + start_time = time() + _loss = function(θ) + if !isnothing(maxtime) && (time() - start_time) > maxtime + error("Optimization halted: time limit exceeded") + end + θ_vec = [θ] + x = cache.f(θ_vec, cache.p) + x = isa(x, Tuple) ? x : (x,) + opt_state = Optimization.OptimizationState(u = θ_vec, objective = x[1]) + if cache.callback(opt_state, x...) + error("Optimization halted by callback") + end + return x[1] + end + kwargs = Dict{Symbol, Any}() + if cache.opt.method == "bounded" + if !isnothing(cache.lb) && !isnothing(cache.ub) + kwargs[:bounds] = (cache.lb[1], cache.ub[1]) + else + throw(ArgumentError("Bounded method requires bounds")) + end + end + _merge_solver_kwargs!(kwargs, cache.solver_args) + t0 = time() + result = nothing + try + result = scipy.optimize.minimize_scalar( + _loss, + method = cache.opt.method; + kwargs... + ) + catch e + if e isa PythonCall.Core.PyException + py_msg = sprint(showerror, e) + if occursin("Optimization halted", py_msg) + throw(ErrorException(py_msg)) + else + throw(ErrorException("SciPy optimization failed: $py_msg")) + end + else + rethrow(e) + end + end + if isnothing(result) + throw(ErrorException("Optimization failed to return a result")) + end + t1 = time() + if pyis(result.x, pybuiltins.None) + minimizer = [NaN] + else + minimizer = [safe_to_float(result.x)] + end + minimum = pyis(result.fun, pybuiltins.None) ? NaN : safe_to_float(result.fun) + py_success = pyconvert(Bool, pybool(result.success)) + if cache.sense === Optimization.MaxSense + minimum = -minimum + end + retcode = py_success ? SciMLBase.ReturnCode.Success : SciMLBase.ReturnCode.Failure + stats = extract_stats(result, t1 - t0) + return SciMLBase.build_solution(cache, cache.opt, minimizer, minimum; + original = result, + retcode = retcode, + stats = stats) +end + +function SciMLBase.__solve(cache::OptimizationCache{F,RC,LB,UB,LC,UC,S,O,D,P,C}) where + {F,RC,LB,UB,LC,UC,S,O<:ScipyLeastSquares,D,P,C} + _residuals = nothing + if hasfield(typeof(cache.f), :f) && (cache.f.f isa ResidualObjective) + real_res = (cache.f.f)::ResidualObjective + _residuals = function(θ) + θ_julia = ensure_julia_array(θ, eltype(cache.u0)) + return real_res.residual(θ_julia, cache.p) + end + else + _residuals = _create_loss(cache; vector_output=true) + end + kwargs = Dict{Symbol, Any}() + kwargs[:method] = cache.opt.method + kwargs[:loss] = cache.opt.loss + if !isnothing(cache.lb) && !isnothing(cache.ub) && cache.opt.method in ["trf", "dogbox"] + kwargs[:bounds] = (cache.lb, cache.ub) + elseif cache.opt.method == "lm" && (!isnothing(cache.lb) || !isnothing(cache.ub)) + @warn "Method 'lm' does not support bounds. Ignoring bounds." + end + kwargs[:jac] = "2-point" + maxiters = Optimization._check_and_convert_maxiters(cache.solver_args.maxiters) + if !isnothing(maxiters) + kwargs[:max_nfev] = maxiters + end + if !isnothing(cache.solver_args.abstol) + kwargs[:ftol] = cache.solver_args.abstol + end + if !isnothing(cache.solver_args.reltol) + kwargs[:gtol] = cache.solver_args.reltol + end + _merge_solver_kwargs!(kwargs, cache.solver_args) + t0 = time() + result = nothing + try + result = scipy.optimize.least_squares( + _residuals, + cache.u0; + kwargs... + ) + catch e + if e isa PythonCall.Core.PyException + py_msg = sprint(showerror, e) + if occursin("Optimization halted by callback", py_msg) + throw(ErrorException("Optimization halted by callback")) + elseif occursin("Optimization halted: time limit exceeded", py_msg) + throw(ErrorException("Optimization halted: time limit exceeded")) + else + throw(ErrorException("SciPy optimization failed: $py_msg")) + end + else + rethrow(e) + end + end + if isnothing(result) + throw(ErrorException("Optimization failed to return a result")) + end + t1 = time() + if pyis(result.x, pybuiltins.None) + minimizer = fill(NaN, length(cache.u0)) + else + minimizer = pyconvert(Vector{eltype(cache.u0)}, result.x) + end + minimum = safe_to_float(result.cost) + py_success = pyconvert(Bool, pybool(result.success)) + py_message = safe_get_message(result) + status = 0 + if pyhasattr(result, "status") + try + status = pyconvert(Int, result.status) + catch + end + end + retcode = scipy_status_to_retcode(status, py_success) + if retcode != SciMLBase.ReturnCode.Success + @debug "ScipyLeastSquares convergence: $(py_message)" + end + stats = extract_stats(result, t1 - t0) + return SciMLBase.build_solution(cache, cache.opt, minimizer, minimum; + original = result, + retcode = retcode, + stats = stats) +end + +function SciMLBase.__solve(cache::OptimizationCache{F,RC,LB,UB,LC,UC,S,O,D,P,C}) where + {F,RC,LB,UB,LC,UC,S,O<:ScipyRootScalar,D,P,C} + x0 = cache.u0[1] + maxtime = get(cache.solver_args, :maxtime, nothing) + start_time = time() + _func = function(θ) + if !isnothing(maxtime) && (time() - start_time) > maxtime + error("Optimization halted: time limit exceeded") + end + θ_vec = [θ] + x = cache.f(θ_vec, cache.p) + x = isa(x, Tuple) ? x : (x,) + opt_state = Optimization.OptimizationState(u = θ_vec, objective = x[1]) + if cache.callback(opt_state, x...) + error("Optimization halted by callback") + end + return x[1] + end + kwargs = Dict{Symbol, Any}() + bracketing_methods = ["brentq", "brenth", "bisect", "ridder"] + is_bracketing = cache.opt.method in bracketing_methods + if is_bracketing + if !isnothing(cache.lb) && !isnothing(cache.ub) + kwargs[:bracket] = pytuple([cache.lb[1], cache.ub[1]]) + else + throw(ArgumentError("Method $(cache.opt.method) requires bracket (bounds)")) + end + else + kwargs[:x0] = x0 + end + if cache.opt.method == "newton" && !isnothing(cache.f.grad) + _fprime = function(θ) + grad = zeros(eltype(cache.u0), 1) + cache.f.grad(grad, [θ], cache.p) + return grad[1] + end + kwargs[:fprime] = _fprime + elseif cache.opt.method == "halley" + if !isnothing(cache.f.grad) && !isnothing(cache.f.hess) + _fprime = function(θ) + grad = zeros(eltype(cache.u0), 1) + cache.f.grad(grad, [θ], cache.p) + return grad[1] + end + _fprime2 = function(θ) + hess = zeros(eltype(cache.u0), 1, 1) + cache.f.hess(hess, [θ], cache.p) + return hess[1, 1] + end + kwargs[:fprime] = _fprime + kwargs[:fprime2] = _fprime2 + else + throw(ArgumentError("Method 'halley' requires both gradient and Hessian")) + end + end + _merge_solver_kwargs!(kwargs, cache.solver_args) + t0 = time() + result = nothing + try + result = scipy.optimize.root_scalar( + _func; + method = cache.opt.method, + kwargs... + ) + catch e + if e isa PythonCall.Core.PyException + py_msg = sprint(showerror, e) + if occursin("Optimization halted", py_msg) + throw(ErrorException(py_msg)) + else + throw(ErrorException("SciPy root finding failed: $py_msg")) + end + else + rethrow(e) + end + end + if isnothing(result) + throw(ErrorException("Root finding failed to return a result")) + end + t1 = time() + if pyis(result.root, pybuiltins.None) + minimizer = [NaN] + root_julia = NaN + minimum = NaN + else + val = safe_to_float(result.root) + minimizer = [val] + root_julia = val + minimum = abs(_func(root_julia)) + end + converged = pyhasattr(result, "converged") ? pyconvert(Bool, pybool(result.converged)) : abs(_func(root_julia)) < 1e-10 + retcode = converged ? SciMLBase.ReturnCode.Success : SciMLBase.ReturnCode.Failure + stats_dict = Dict{Symbol, Any}(:time => t1 - t0) + if pyhasattr(result, "iterations") + try stats_dict[:iterations] = pyconvert(Int, result.iterations) catch; end + end + if pyhasattr(result, "function_calls") + try stats_dict[:fevals] = pyconvert(Int, result.function_calls) catch; end + end + stats = Optimization.OptimizationStats(; stats_dict...) + return SciMLBase.build_solution(cache, cache.opt, minimizer, minimum; + original = result, + retcode = retcode, + stats = stats) +end + +function SciMLBase.__solve(cache::OptimizationCache{F,RC,LB,UB,LC,UC,S,O,D,P,C}) where + {F,RC,LB,UB,LC,UC,S,O<:ScipyRoot,D,P,C} + _func = _create_loss(cache, vector_output=true) + kwargs = Dict{Symbol, Any}() + kwargs[:method] = cache.opt.method + if !isnothing(cache.f.grad) && cache.opt.method in ["hybr", "lm"] + _jac = function(θ) + θ_julia = ensure_julia_array(θ, eltype(cache.u0)) + fval = cache.f(θ_julia, cache.p) + if isa(fval, Tuple) + fval = fval[1] + end + if isa(fval, Number) + fval = [fval] + end + m = length(fval) + n = length(θ_julia) + jac = zeros(eltype(cache.u0), m, n) + cache.f.grad(jac, θ_julia, cache.p) + return jac + end + kwargs[:jac] = _jac + end + if isa(cache.solver_args, NamedTuple) + _merge_solver_kwargs!(kwargs, cache.solver_args) + end + t0 = time() + result = nothing + try + result = scipy.optimize.root( + _func, + cache.u0; + kwargs... + ) + catch e + if e isa PythonCall.Core.PyException + py_msg = sprint(showerror, e) + if occursin("Optimization halted by callback", py_msg) + throw(ErrorException("Optimization halted by callback")) + elseif occursin("Optimization halted: time limit exceeded", py_msg) + throw(ErrorException("Optimization halted: time limit exceeded")) + else + throw(ErrorException("SciPy root finding failed: $py_msg")) + end + else + rethrow(e) + end + end + if isnothing(result) + throw(ErrorException("Root finding failed to return a result")) + end + t1 = time() + if pyis(result.x, pybuiltins.None) + minimizer = fill(NaN, length(cache.u0)) + else + minimizer = pyconvert(Vector{eltype(cache.u0)}, result.x) + end + fun_val = pyconvert(Vector{Float64}, result.fun) + minimum = sum(abs2, fun_val) + py_success = pyconvert(Bool, pybool(result.success)) + py_message = safe_get_message(result) + retcode = py_success ? SciMLBase.ReturnCode.Success : SciMLBase.ReturnCode.Failure + if retcode != SciMLBase.ReturnCode.Success + @debug "ScipyRoot convergence: $(py_message)" + end + stats = extract_stats(result, t1 - t0) + return SciMLBase.build_solution(cache, cache.opt, minimizer, minimum; + original = result, + retcode = retcode, + stats = stats) +end + +function SciMLBase.__solve(cache::OptimizationCache{F,RC,LB,UB,LC,UC,S,O,D,P,C}) where + {F,RC,LB,UB,LC,UC,S,O<:ScipyLinprog,D,P,C} + c = cache.f(cache.u0, cache.p) + if isa(c, Tuple) + c = c[1] + end + if isa(c, Number) + c = [c] + end + bounds = nothing + if !isnothing(cache.lb) || !isnothing(cache.ub) + n = length(cache.u0) + lb = isnothing(cache.lb) ? fill(-Inf, n) : cache.lb + ub = isnothing(cache.ub) ? fill(Inf, n) : cache.ub + if length(lb) != n + lb_new = fill(-Inf, n) + lb_new[1:min(length(lb), n)] .= lb[1:min(length(lb), n)] + lb = lb_new + end + if length(ub) != n + ub_new = fill(Inf, n) + ub_new[1:min(length(ub), n)] .= ub[1:min(length(ub), n)] + ub = ub_new + end + bounds_list = [] + for i in 1:n + lb_val = isfinite(lb[i]) ? lb[i] : nothing + ub_val = isfinite(ub[i]) ? ub[i] : nothing + push!(bounds_list, (lb_val, ub_val)) + end + bounds = pylist(bounds_list) + end + # Allow users to pass constraint matrices via solver kwargs + A_ub = get(cache.solver_args, :A_ub, nothing) + b_ub = get(cache.solver_args, :b_ub, nothing) + A_eq = get(cache.solver_args, :A_eq, nothing) + b_eq = get(cache.solver_args, :b_eq, nothing) + if !(isnothing(A_ub) == isnothing(b_ub)) + throw(ArgumentError("Both A_ub and b_ub must be provided together")) + end + if !(isnothing(A_eq) == isnothing(b_eq)) + throw(ArgumentError("Both A_eq and b_eq must be provided together")) + end + maxiters = Optimization._check_and_convert_maxiters(cache.solver_args.maxiters) + options = nothing + if !isnothing(maxiters) + options = pydict(Dict("maxiter" => maxiters)) + end + t0 = time() + result = nothing + try + result = scipy.optimize.linprog( + c, + A_ub = A_ub, + b_ub = b_ub, + A_eq = A_eq, + b_eq = b_eq, + bounds = bounds, + method = cache.opt.method, + options = options + ) + catch e + if e isa PythonCall.Core.PyException + py_msg = sprint(showerror, e) + throw(ErrorException("SciPy linear programming failed: $py_msg")) + else + rethrow(e) + end + end + if isnothing(result) + throw(ErrorException("Linear programming failed to return a result")) + end + t1 = time() + if pyis(result.x, pybuiltins.None) + minimizer = fill(NaN, length(cache.u0)) + else + minimizer = pyconvert(Vector{eltype(cache.u0)}, result.x) + end + minimum = pyis(result.fun, pybuiltins.None) ? NaN : safe_to_float(result.fun) + py_success = pyconvert(Bool, pybool(result.success)) + py_message = safe_get_message(result) + if cache.sense === Optimization.MaxSense + minimum = -minimum + end + status = 0 + if pyhasattr(result, "status") + try + status = pyconvert(Int, result.status) + catch + end + end + retcode = scipy_status_to_retcode(status, py_success) + if retcode != SciMLBase.ReturnCode.Success + @debug "ScipyLinprog convergence: $(py_message)" + end + stats = extract_stats(result, t1 - t0) + return SciMLBase.build_solution(cache, cache.opt, minimizer, minimum; + original = result, + retcode = retcode, + stats = stats) +end + +function SciMLBase.__solve(cache::OptimizationCache{F,RC,LB,UB,LC,UC,S,O,D,P,C}) where + {F,RC,LB,UB,LC,UC,S,O<:ScipyMilp,D,P,C} + c = cache.f(cache.u0, cache.p) + if isa(c, Tuple) + c = c[1] + end + if isa(c, Number) + c = [c] + end + n = length(c) + lb = isnothing(cache.lb) ? fill(-Inf, n) : copy(cache.lb) + ub = isnothing(cache.ub) ? fill(Inf, n) : copy(cache.ub) + if length(lb) != n + lb_new = fill(-Inf, n) + lb_new[1:min(length(lb), n)] .= lb[1:min(length(lb), n)] + lb = lb_new + end + if length(ub) != n + ub_new = fill(Inf, n) + ub_new[1:min(length(ub), n)] .= ub[1:min(length(ub), n)] + ub = ub_new + end + bounds = scipy.optimize.Bounds(lb, ub) + integrality = get(cache.solver_args, :integrality, nothing) + A = get(cache.solver_args, :A, nothing) + lb_con = get(cache.solver_args, :lb_con, nothing) + ub_con = get(cache.solver_args, :ub_con, nothing) + constraints = nothing + if !(isnothing(A) && isnothing(lb_con) && isnothing(ub_con)) + if any(isnothing.((A, lb_con, ub_con))) + throw(ArgumentError("A, lb_con, and ub_con must all be provided for linear constraints")) + end + keep_feasible_flag = get(cache.solver_args, :keep_feasible, false) + constraints = scipy.optimize.LinearConstraint(A, lb_con, ub_con, keep_feasible = keep_feasible_flag) + end + t0 = time() + result = nothing + try + result = scipy.optimize.milp( + c = c, + integrality = integrality, + bounds = bounds, + constraints = constraints, + options = nothing + ) + catch e + if e isa PythonCall.Core.PyException + py_msg = sprint(showerror, e) + throw(ErrorException("SciPy MILP failed: $py_msg")) + else + rethrow(e) + end + end + if isnothing(result) + throw(ErrorException("MILP failed to return a result")) + end + t1 = time() + if pyis(result.x, pybuiltins.None) + minimizer = fill(NaN, length(cache.u0)) + else + minimizer = pyconvert(Vector{eltype(cache.u0)}, result.x) + end + minimum = pyis(result.fun, pybuiltins.None) ? NaN : safe_to_float(result.fun) + py_success = pyconvert(Bool, pybool(result.success)) + py_message = safe_get_message(result) + if cache.sense === Optimization.MaxSense + minimum = -minimum + end + retcode = py_success ? SciMLBase.ReturnCode.Success : SciMLBase.ReturnCode.Failure + if retcode != SciMLBase.ReturnCode.Success + @debug "ScipyMilp convergence: $(py_message)" + end + stats = extract_stats(result, t1 - t0) + return SciMLBase.build_solution(cache, cache.opt, minimizer, minimum; + original = result, + retcode = retcode, + stats = stats) +end + +function SciMLBase.__solve(cache::OptimizationCache{F,RC,LB,UB,LC,UC,S,O,D,P,C}) where + {F,RC,LB,UB,LC,UC,S,O<:ScipyDifferentialEvolution,D,P,C} + _loss = _create_loss(cache) + bounds = _build_bounds(cache.lb, cache.ub) + maxiters = Optimization._check_and_convert_maxiters(cache.solver_args.maxiters) + de_kwargs = Dict{Symbol, Any}() + de_kwargs[:maxiter] = isnothing(maxiters) ? 1000 : maxiters + de_kwargs[:popsize] = 15 + de_kwargs[:atol] = 0.0 + de_kwargs[:tol] = 0.01 + de_kwargs[:mutation] = (0.5, 1) + de_kwargs[:recombination] = 0.7 + de_kwargs[:polish] = true + de_kwargs[:init] = "latinhypercube" + de_kwargs[:updating] = "immediate" + de_kwargs[:workers] = 1 + _merge_solver_kwargs!(de_kwargs, cache.solver_args) + t0 = time() + result = nothing + try + result = scipy.optimize.differential_evolution( + _loss, + bounds; + de_kwargs... + ) + catch e + if e isa PythonCall.Core.PyException + py_msg = sprint(showerror, e) + if occursin("Optimization halted by callback", py_msg) + throw(ErrorException("Optimization halted by callback")) + elseif occursin("Optimization halted: time limit exceeded", py_msg) + throw(ErrorException("Optimization halted: time limit exceeded")) + else + throw(ErrorException("SciPy optimization failed: $py_msg")) + end + else + rethrow(e) + end + end + if isnothing(result) + throw(ErrorException("Optimization failed to return a result")) + end + t1 = time() + if pyis(result.x, pybuiltins.None) + minimizer = fill(NaN, length(cache.u0)) + else + minimizer = pyconvert(Vector{eltype(cache.u0)}, result.x) + end + minimum = safe_to_float(result.fun) + py_success = pyconvert(Bool, pybool(result.success)) + py_message = safe_get_message(result) + if cache.sense === Optimization.MaxSense + minimum = -minimum + end + retcode = py_success ? SciMLBase.ReturnCode.Success : SciMLBase.ReturnCode.Failure + if retcode != SciMLBase.ReturnCode.Success + @debug "ScipyDifferentialEvolution convergence: $(py_message)" + end + stats = extract_stats(result, t1 - t0) + return SciMLBase.build_solution(cache, cache.opt, minimizer, minimum; + original = result, + retcode = retcode, + stats = stats) +end + +function SciMLBase.__solve(cache::OptimizationCache{F,RC,LB,UB,LC,UC,S,O,D,P,C}) where + {F,RC,LB,UB,LC,UC,S,O<:ScipyBasinhopping,D,P,C} + _loss = _create_loss(cache) + maxiters = Optimization._check_and_convert_maxiters(cache.solver_args.maxiters) + bh_kwargs = Dict{Symbol, Any}() + bh_kwargs[:niter] = isnothing(maxiters) ? 100 : maxiters + bh_kwargs[:T] = 1.0 + bh_kwargs[:stepsize] = 0.5 + bh_kwargs[:interval] = 50 + _merge_solver_kwargs!(bh_kwargs, cache.solver_args) + t0 = time() + result = nothing + try + result = scipy.optimize.basinhopping( + _loss, + cache.u0; + bh_kwargs... + ) + catch e + if e isa PythonCall.Core.PyException + py_msg = sprint(showerror, e) + if occursin("Optimization halted by callback", py_msg) + throw(ErrorException("Optimization halted by callback")) + elseif occursin("Optimization halted: time limit exceeded", py_msg) + throw(ErrorException("Optimization halted: time limit exceeded")) + else + throw(ErrorException("SciPy optimization failed: $py_msg")) + end + else + rethrow(e) + end + end + if isnothing(result) + throw(ErrorException("Optimization failed to return a result")) + end + t1 = time() + if pyis(result.x, pybuiltins.None) + minimizer = fill(NaN, length(cache.u0)) + else + minimizer = pyconvert(Vector{eltype(cache.u0)}, result.x) + end + minimum = safe_to_float(result.fun) + lowest_result = result.lowest_optimization_result + py_success = pyconvert(Bool, pybool(lowest_result.success)) + py_message = safe_get_message(lowest_result) + if cache.sense === Optimization.MaxSense + minimum = -minimum + end + retcode = py_success ? SciMLBase.ReturnCode.Success : SciMLBase.ReturnCode.Failure + if retcode != SciMLBase.ReturnCode.Success + @debug "ScipyBasinhopping convergence: $(py_message)" + end + stats = extract_stats(lowest_result, t1 - t0) + return SciMLBase.build_solution(cache, cache.opt, minimizer, minimum; + original = result, + retcode = retcode, + stats = stats) +end + +function SciMLBase.__solve(cache::OptimizationCache{F,RC,LB,UB,LC,UC,S,O,D,P,C}) where + {F,RC,LB,UB,LC,UC,S,O<:ScipyDualAnnealing,D,P,C} + _loss = _create_loss(cache) + bounds = _build_bounds(cache.lb, cache.ub) + da_kwargs = Dict{Symbol, Any}() + da_kwargs[:maxiter] = begin + mi = Optimization._check_and_convert_maxiters(cache.solver_args.maxiters) + isnothing(mi) ? 1000 : mi + end + da_kwargs[:initial_temp] = 5230.0 + da_kwargs[:restart_temp_ratio] = 2e-5 + da_kwargs[:visit] = 2.62 + da_kwargs[:accept] = -5.0 + da_kwargs[:maxfun] = 1e7 + da_kwargs[:no_local_search] = false + _merge_solver_kwargs!(da_kwargs, cache.solver_args) + t0 = time() + result = nothing + try + result = scipy.optimize.dual_annealing( + _loss, + bounds; + da_kwargs... + ) + catch e + if e isa PythonCall.Core.PyException + py_msg = sprint(showerror, e) + if occursin("Optimization halted by callback", py_msg) + throw(ErrorException("Optimization halted by callback")) + elseif occursin("Optimization halted: time limit exceeded", py_msg) + throw(ErrorException("Optimization halted: time limit exceeded")) + else + throw(ErrorException("SciPy optimization failed: $py_msg")) + end + else + rethrow(e) + end + end + if isnothing(result) + throw(ErrorException("Optimization failed to return a result")) + end + t1 = time() + if pyis(result.x, pybuiltins.None) + minimizer = fill(NaN, length(cache.u0)) + else + minimizer = pyconvert(Vector{eltype(cache.u0)}, result.x) + end + minimum = safe_to_float(result.fun) + py_success = pyconvert(Bool, pybool(result.success)) + py_message = safe_get_message(result) + if cache.sense === Optimization.MaxSense + minimum = -minimum + end + retcode = py_success ? SciMLBase.ReturnCode.Success : SciMLBase.ReturnCode.Failure + if retcode != SciMLBase.ReturnCode.Success + @debug "ScipyDualAnnealing convergence: $(py_message)" + end + stats = extract_stats(result, t1 - t0) + return SciMLBase.build_solution(cache, cache.opt, minimizer, minimum; + original = result, + retcode = retcode, + stats = stats) +end + +function SciMLBase.__solve(cache::OptimizationCache{F,RC,LB,UB,LC,UC,S,O,D,P,C}) where + {F,RC,LB,UB,LC,UC,S,O<:ScipyShgo,D,P,C} + local cons_cache = nothing + if !isnothing(cache.f.cons) && !isnothing(cache.lcons) + cons_cache = zeros(eltype(cache.u0), length(cache.lcons)) + end + _loss = _create_loss(cache) + bounds = _build_bounds(cache.lb, cache.ub) + constraints = nothing + if !isnothing(cons_cache) + cons_list = [] + _cons_func = function(θ) + θ_julia = ensure_julia_array(θ, eltype(cache.u0)) + cons_cache .= zero(eltype(cons_cache)) + if hasmethod(cache.f.cons, Tuple{typeof(cons_cache), typeof(θ_julia), typeof(cache.p)}) + cache.f.cons(cons_cache, θ_julia, cache.p) + else + cache.f.cons(cons_cache, θ_julia) + end + return cons_cache + end + for i in 1:length(cache.lcons) + if isfinite(cache.lcons[i]) + cons_func_i = let i=i, _cons_func=_cons_func + θ -> _cons_func(θ)[i] - cache.lcons[i] + end + push!(cons_list, pydict(Dict("type" => "ineq", "fun" => cons_func_i))) + end + end + for i in 1:length(cache.ucons) + if isfinite(cache.ucons[i]) + cons_func_i = let i=i, _cons_func=_cons_func + θ -> cache.ucons[i] - _cons_func(θ)[i] + end + push!(cons_list, pydict(Dict("type" => "ineq", "fun" => cons_func_i))) + end + end + constraints = pylist(cons_list) + end + shgo_kwargs = Dict{Symbol, Any}() + shgo_kwargs[:n] = 100 + shgo_kwargs[:iters] = 1 + shgo_kwargs[:sampling_method] = "simplicial" + _merge_solver_kwargs!(shgo_kwargs, cache.solver_args) + t0 = time() + result = nothing + try + result = scipy.optimize.shgo( + _loss, + bounds; + args = (), + constraints = constraints, + shgo_kwargs... + ) + catch e + if e isa PythonCall.Core.PyException + py_msg = sprint(showerror, e) + if occursin("Optimization halted by callback", py_msg) + throw(ErrorException("Optimization halted by callback")) + elseif occursin("Optimization halted: time limit exceeded", py_msg) + throw(ErrorException("Optimization halted: time limit exceeded")) + else + throw(ErrorException("SciPy optimization failed: $py_msg")) + end + else + rethrow(e) + end + end + if isnothing(result) + throw(ErrorException("Optimization failed to return a result")) + end + t1 = time() + if pyis(result.x, pybuiltins.None) + minimizer = fill(NaN, length(cache.u0)) + else + minimizer = pyconvert(Vector{eltype(cache.u0)}, result.x) + end + minimum = safe_to_float(result.fun) + py_success = pyconvert(Bool, pybool(result.success)) + py_message = safe_get_message(result) + if cache.sense === Optimization.MaxSense + minimum = -minimum + end + retcode = py_success ? SciMLBase.ReturnCode.Success : SciMLBase.ReturnCode.Failure + if retcode != SciMLBase.ReturnCode.Success + @debug "ScipyShgo convergence: $(py_message)" + end + stats = extract_stats(result, t1 - t0) + return SciMLBase.build_solution(cache, cache.opt, minimizer, minimum; + original = result, + retcode = retcode, + stats = stats) +end + +function SciMLBase.__solve(cache::OptimizationCache{F,RC,LB,UB,LC,UC,S,O,D,P,C}) where + {F,RC,LB,UB,LC,UC,S,O<:ScipyDirect,D,P,C} + _loss = _create_loss(cache) + bounds = _build_bounds(cache.lb, cache.ub) + maxiters = Optimization._check_and_convert_maxiters(cache.solver_args.maxiters) + direct_kwargs = Dict{Symbol, Any}() + direct_kwargs[:eps] = 0.0001 + direct_kwargs[:maxiter] = isnothing(maxiters) ? 1000 : maxiters + direct_kwargs[:locally_biased] = true + direct_kwargs[:vol_tol] = 1e-16 + direct_kwargs[:len_tol] = 1e-6 + _merge_solver_kwargs!(direct_kwargs, cache.solver_args) + t0 = time() + result = nothing + try + result = scipy.optimize.direct( + _loss, + bounds; + direct_kwargs... + ) + catch e + if e isa PythonCall.Core.PyException + py_msg = sprint(showerror, e) + if occursin("Optimization halted by callback", py_msg) + throw(ErrorException("Optimization halted by callback")) + elseif occursin("Optimization halted: time limit exceeded", py_msg) + throw(ErrorException("Optimization halted: time limit exceeded")) + else + throw(ErrorException("SciPy optimization failed: $py_msg")) + end + else + rethrow(e) + end + end + if isnothing(result) + throw(ErrorException("Optimization failed to return a result")) + end + t1 = time() + if pyis(result.x, pybuiltins.None) + minimizer = fill(NaN, length(cache.u0)) + else + minimizer = pyconvert(Vector{eltype(cache.u0)}, result.x) + end + minimum = safe_to_float(result.fun) + py_success = pyconvert(Bool, pybool(result.success)) + py_message = safe_get_message(result) + if cache.sense === Optimization.MaxSense + minimum = -minimum + end + retcode = py_success ? SciMLBase.ReturnCode.Success : SciMLBase.ReturnCode.Failure + if retcode != SciMLBase.ReturnCode.Success + @debug "ScipyDirect convergence: $(py_message)" + end + stats = extract_stats(result, t1 - t0) + return SciMLBase.build_solution(cache, cache.opt, minimizer, minimum; + original = result, + retcode = retcode, + stats = stats) +end + +function SciMLBase.__solve(cache::OptimizationCache{F,RC,LB,UB,LC,UC,S,O,D,P,C}) where + {F,RC,LB,UB,LC,UC,S,O<:ScipyBrute,D,P,C} + _loss = _create_loss(cache) + ranges = _build_bounds(cache.lb, cache.ub) + brute_kwargs = Dict{Symbol, Any}() + brute_kwargs[:Ns] = 20 + brute_kwargs[:full_output] = true + brute_kwargs[:finish] = scipy.optimize.fmin + brute_kwargs[:workers] = 1 + _merge_solver_kwargs!(brute_kwargs, cache.solver_args) + t0 = time() + result = nothing + try + result = scipy.optimize.brute( + _loss, + ranges; + brute_kwargs... + ) + catch e + if e isa PythonCall.Core.PyException + py_msg = sprint(showerror, e) + if occursin("Optimization halted by callback", py_msg) + throw(ErrorException("Optimization halted by callback")) + elseif occursin("Optimization halted: time limit exceeded", py_msg) + throw(ErrorException("Optimization halted: time limit exceeded")) + else + throw(ErrorException("SciPy optimization failed: $py_msg")) + end + else + rethrow(e) + end + end + if isnothing(result) + throw(ErrorException("Optimization failed to return a result")) + end + t1 = time() + if pyis(result[0], pybuiltins.None) + minimizer = fill(NaN, length(cache.u0)) + else + minimizer = pyconvert(Vector{eltype(cache.u0)}, result[0]) + end + minimum = safe_to_float(result[1]) + if cache.sense === Optimization.MaxSense + minimum = -minimum + end + retcode = SciMLBase.ReturnCode.Success + stats = Optimization.OptimizationStats(; time = t1 - t0) + return SciMLBase.build_solution(cache, cache.opt, minimizer, minimum; + original = result, + retcode = retcode, + stats = stats) +end + +export ScipyMinimize, ScipyNelderMead, ScipyPowell, ScipyCG, ScipyBFGS, ScipyNewtonCG, + ScipyLBFGSB, ScipyTNC, ScipyCOBYLA, ScipyCOBYQA, ScipySLSQP, ScipyTrustConstr, + ScipyDogleg, ScipyTrustNCG, ScipyTrustKrylov, ScipyTrustExact, + ScipyMinimizeScalar, ScipyBrent, ScipyBounded, ScipyGolden, + ScipyLeastSquares, ScipyLeastSquaresTRF, ScipyLeastSquaresDogbox, ScipyLeastSquaresLM, + ScipyRootScalar, ScipyRoot, ScipyLinprog, ScipyMilp, + ScipyDifferentialEvolution, ScipyBasinhopping, ScipyDualAnnealing, + ScipyShgo, ScipyDirect, ScipyBrute + +# Wrap the user's Julia objective so it matches what SciPy expects. +function _create_loss(cache; vector_output::Bool = false) + maxtime = get(cache.solver_args, :maxtime, nothing) + start_time = !isnothing(maxtime) ? time() : 0.0 + if vector_output + return function (θ) + if !isnothing(maxtime) && (time() - start_time) > maxtime + error("Optimization halted: time limit exceeded") + end + θ_julia = ensure_julia_array(θ, eltype(cache.u0)) + x = cache.f(θ_julia, cache.p) + if isa(x, Tuple) + x = x + elseif isa(x, Number) + x = (x,) + end + opt_state = Optimization.OptimizationState(u = θ_julia, objective = sum(abs2, x)) + if cache.callback(opt_state, x...) + error("Optimization halted by callback") + end + + arr = cache.sense === Optimization.MaxSense ? -x : x + return arr + end + else + return function (θ) + if !isnothing(maxtime) && (time() - start_time) > maxtime + error("Optimization halted: time limit exceeded") + end + θ_julia = ensure_julia_array(θ, eltype(cache.u0)) + x = cache.f(θ_julia, cache.p) + if isa(x, Tuple) + x = x + elseif isa(x, Number) + x = (x,) + end + opt_state = Optimization.OptimizationState(u = θ_julia, objective = x[1]) + if cache.callback(opt_state, x...) + error("Optimization halted by callback") + end + return cache.sense === Optimization.MaxSense ? -x[1] : x[1] + end + end +end + +# These solver-args are handled specially elsewhere, so we skip them here. +const _DEFAULT_EXCLUDE = ( + :maxiters, :maxtime, :abstol, :reltol, :callback, :progress, :cons_tol, + :jac_sparsity, :keep_feasible, :hess_update +) + +# Moving the remaining kwargs into a Dict that we pass straight to SciPy. +function _merge_solver_kwargs!(dest::AbstractDict, solver_args; exclude = _DEFAULT_EXCLUDE) + if isa(solver_args, NamedTuple) + for (k, v) in pairs(solver_args) + k in exclude && continue + isnothing(v) && continue + dest[convert(keytype(dest), k)] = v + end + end + return dest +end + +function _build_bounds(lb::AbstractVector, ub::AbstractVector) + return pylist([pytuple([lb[i], ub[i]]) for i in eachindex(lb)]) +end + +struct ResidualObjective{R} + residual::R +end + +(r::ResidualObjective)(u, p) = sum(abs2, r.residual(u, p)) + +end + diff --git a/lib/OptimizationSciPy/test/runtests.jl b/lib/OptimizationSciPy/test/runtests.jl new file mode 100644 index 000000000..c178110fc --- /dev/null +++ b/lib/OptimizationSciPy/test/runtests.jl @@ -0,0 +1,486 @@ +using OptimizationSciPy, Optimization, Zygote, ReverseDiff, ForwardDiff +using Test, Random +using Optimization.SciMLBase: ReturnCode, NonlinearLeastSquaresProblem +using PythonCall + +function rosenbrock(x, p) + (p[1] - x[1])^2 + p[2] * (x[2] - x[1]^2)^2 +end + +function rosenbrock_hess(H, x, p) + H[1,1] = 2 - 400*p[2]*x[2] + 1200*p[2]*x[1]^2 + H[1,2] = -400*p[2]*x[1] + H[2,1] = -400*p[2]*x[1] + H[2,2] = 200*p[2] + return nothing +end + +@testset "OptimizationSciPy.jl" begin + x0 = zeros(2) + _p = [1.0, 100.0] + l1 = rosenbrock(x0, _p) + + @testset "MaxSense" begin + optprob = OptimizationFunction((x, p) -> -rosenbrock(x, p), Optimization.AutoZygote()) + prob = OptimizationProblem(optprob, x0, _p; sense = Optimization.MaxSense) + sol = solve(prob, ScipyNelderMead()) + @test sol.retcode == ReturnCode.Success + @test 10 * sol.objective < l1 + end + + @testset "unconstrained with gradient" begin + optprob = OptimizationFunction(rosenbrock, Optimization.AutoZygote()) + prob = OptimizationProblem(optprob, x0, _p) + sol = solve(prob, ScipyBFGS()) + @test sol.retcode == ReturnCode.Success + @test 10 * sol.objective < l1 + sol = solve(prob, ScipyLBFGSB()) + @test sol.retcode == ReturnCode.Success + @test 10 * sol.objective < l1 + end + + @testset "bounded" begin + optprob = OptimizationFunction(rosenbrock, Optimization.AutoZygote()) + prob = OptimizationProblem(optprob, x0, _p, lb = [-1.0, -1.0], ub = [0.8, 0.8]) + sol = solve(prob, ScipyLBFGSB()) + @test sol.retcode == ReturnCode.Success + @test 10 * sol.objective < l1 + end + + @testset "global optimization" begin + optprob = OptimizationFunction(rosenbrock, Optimization.AutoZygote()) + prob = OptimizationProblem(optprob, x0, _p, lb = [-1.0, -1.0], ub = [0.8, 0.8]) + sol = solve(prob, ScipyDifferentialEvolution(), maxiters = 100) + @test sol.retcode == ReturnCode.Success + @test 10 * sol.objective < l1 + sol = solve(prob, ScipyBasinhopping(), maxiters = 50) + @test sol.retcode == ReturnCode.Success + @test 10 * sol.objective < l1 + sol = solve(prob, ScipyDualAnnealing(), maxiters = 100) + @test sol.retcode == ReturnCode.Success + @test 10 * sol.objective < l1 + sol = solve(prob, ScipyShgo()) + @test sol.retcode == ReturnCode.Success + @test 10 * sol.objective < l1 + sol = solve(prob, ScipyDirect(), maxiters = 1000) + @test sol.retcode in (ReturnCode.Success, ReturnCode.Failure) + if sol.retcode == ReturnCode.Success + @test 10 * sol.objective < l1 + end + sol = solve(prob, ScipyBrute(), Ns = 10) + @test sol.retcode == ReturnCode.Success + @test 10 * sol.objective < l1 + end + + @testset "various methods" begin + optprob = OptimizationFunction(rosenbrock, Optimization.AutoZygote()) + prob = OptimizationProblem(optprob, x0, _p) + sol = solve(prob, ScipyNelderMead()) + @test sol.retcode == ReturnCode.Success + @test 10 * sol.objective < l1 + sol = solve(prob, ScipyPowell()) + @test sol.retcode == ReturnCode.Success + @test 10 * sol.objective < l1 + sol = solve(prob, ScipyCG()) + @test sol.retcode == ReturnCode.Success + @test 10 * sol.objective < l1 + sol = solve(prob, ScipyTNC()) + @test sol.retcode == ReturnCode.Success + @test 10 * sol.objective < l1 + end + + @testset "with Hessian" begin + optf = OptimizationFunction(rosenbrock, Optimization.AutoForwardDiff(); hess = rosenbrock_hess) + prob = OptimizationProblem(optf, x0, _p) + sol = solve(prob, ScipyNewtonCG(), maxiters = 200) + @test sol.retcode == ReturnCode.Success + @test 10 * sol.objective < l1 + end + + @testset "bounded optimization" begin + optprob = OptimizationFunction(rosenbrock, Optimization.AutoZygote()) + prob = OptimizationProblem(optprob, x0, _p, lb = [-1.0, -1.0], ub = [0.8, 0.8]) + sol = solve(prob, ScipyLBFGSB()) + @test sol.retcode == ReturnCode.Success + @test 10 * sol.objective < l1 + sol = solve(prob, ScipyTNC()) + @test sol.retcode == ReturnCode.Success + @test 10 * sol.objective < l1 + end + + @testset "trust region with Hessian" begin + optf_hess = OptimizationFunction(rosenbrock, Optimization.AutoForwardDiff(); hess = rosenbrock_hess) + x0_trust = [0.5, 0.5] + prob = OptimizationProblem(optf_hess, x0_trust, _p) + for method in [ScipyDogleg(), ScipyTrustNCG(), ScipyTrustKrylov(), ScipyTrustExact()] + sol = solve(prob, method, maxiters = 2000) + @test sol.retcode in (ReturnCode.Success, ReturnCode.MaxIters, ReturnCode.Unstable, ReturnCode.Infeasible) + if sol.retcode == ReturnCode.Success + @test 10 * sol.objective < sol.original.fun + end + end + end + + @testset "COBYQA method" begin + optf_no_grad = OptimizationFunction(rosenbrock) + prob = OptimizationProblem(optf_no_grad, x0, _p) + sol = solve(prob, ScipyCOBYQA(), maxiters = 10000) + @test sol.retcode == ReturnCode.Success + @test 10 * sol.objective < l1 + prob_bounded = OptimizationProblem(optf_no_grad, x0, _p, lb = [-1.0, -1.0], ub = [0.8, 0.8]) + sol = solve(prob_bounded, ScipyCOBYQA()) + @test sol.retcode == ReturnCode.Success + cons = (res, x, p) -> res .= [x[1]^2 + x[2]^2 - 1.0] + optf_cons = OptimizationFunction(rosenbrock; cons = cons) + prob_cons = OptimizationProblem(optf_cons, [0.5, 0.5], _p, lcons = [-0.01], ucons = [0.01]) + sol = solve(prob_cons, ScipyCOBYQA()) + @test sol.retcode == ReturnCode.Success + end + + @testset "ScipyMinimizeScalar" begin + f_scalar(x, p) = (x[1] - p[1])^2 + sin(x[1]) + x0_scalar = [0.0] + p_scalar = [2.0] + optf = OptimizationFunction(f_scalar) + prob = OptimizationProblem(optf, x0_scalar, p_scalar) + sol = solve(prob, ScipyBrent()) + @test sol.retcode == ReturnCode.Success + @test length(sol.u) == 1 + @test abs(2*(sol.u[1] - p_scalar[1]) + cos(sol.u[1])) < 1e-6 + sol = solve(prob, ScipyGolden()) + @test sol.retcode == ReturnCode.Success + @test abs(2*(sol.u[1] - p_scalar[1]) + cos(sol.u[1])) < 1e-6 + prob_bounded = OptimizationProblem(optf, x0_scalar, p_scalar, lb = [0.0], ub = [3.0]) + sol = solve(prob_bounded, ScipyBounded()) + @test sol.retcode == ReturnCode.Success + @test 0.0 <= sol.u[1] <= 3.0 + prob_multidim = OptimizationProblem(rosenbrock, x0, _p) + @test_throws ArgumentError solve(prob_multidim, ScipyMinimizeScalar("brent")) + @test_throws ArgumentError solve(prob, ScipyBounded()) + optf_grad = OptimizationFunction(f_scalar, Optimization.AutoZygote()) + prob_grad = OptimizationProblem(optf_grad, x0_scalar, p_scalar) + sol = solve(prob_grad, ScipyBrent()) + @test sol.retcode == ReturnCode.Success + end + + @testset "ScipyRootScalar" begin + f_root(x, p) = x[1]^3 - 2*x[1] - 5 + x0_root = [2.0] + optf = OptimizationFunction(f_root) + prob_bracket = OptimizationProblem(optf, x0_root, nothing, lb = [2.0], ub = [3.0]) + sol = solve(prob_bracket, ScipyRootScalar("brentq")) + @test sol.retcode == ReturnCode.Success + @test abs(f_root(sol.u, nothing)) < 1e-10 + sol = solve(prob_bracket, ScipyRootScalar("brenth")) + @test sol.retcode == ReturnCode.Success + @test abs(f_root(sol.u, nothing)) < 1e-10 + sol = solve(prob_bracket, ScipyRootScalar("bisect")) + @test sol.retcode == ReturnCode.Success + @test abs(f_root(sol.u, nothing)) < 1e-10 + sol = solve(prob_bracket, ScipyRootScalar("ridder")) + @test sol.retcode == ReturnCode.Success + @test abs(f_root(sol.u, nothing)) < 1e-10 + prob_no_bracket = OptimizationProblem(optf, x0_root) + sol = solve(prob_no_bracket, ScipyRootScalar("secant")) + @test sol.retcode == ReturnCode.Success + @test abs(f_root(sol.u, nothing)) < 1e-10 + f_root_grad(g, x, p) = g[1] = 3*x[1]^2 - 2 + optf_grad = OptimizationFunction(f_root; grad = f_root_grad) + prob_newton = OptimizationProblem(optf_grad, x0_root) + sol = solve(prob_newton, ScipyRootScalar("newton")) + @test sol.retcode == ReturnCode.Success + @test abs(f_root(sol.u, nothing)) < 1e-10 + f_root_hess(H, x, p) = H[1,1] = 6*x[1] + optf_halley = OptimizationFunction(f_root; grad = f_root_grad, hess = f_root_hess) + prob_halley = OptimizationProblem(optf_halley, x0_root) + sol = solve(prob_halley, ScipyRootScalar("halley")) + @test sol.retcode == ReturnCode.Success + @test abs(f_root(sol.u, nothing)) < 1e-10 + prob_multidim = OptimizationProblem(rosenbrock, x0, _p) + @test_throws ArgumentError solve(prob_multidim, ScipyRootScalar("brentq")) + @test_throws ArgumentError solve(prob_no_bracket, ScipyRootScalar("brentq")) + end + + @testset "ScipyRoot" begin + function system(x, p) + return [x[1]^2 + x[2]^2 - 1.0, x[2] - x[1]^2] + end + x0_system = [0.5, 0.5] + optf = OptimizationFunction(system) + prob = OptimizationProblem(optf, x0_system) + sol = solve(prob, ScipyRoot("hybr")) + @test sol.retcode == ReturnCode.Success + res = system(sol.u, nothing) + @test all(abs.(res) .< 1e-10) + sol = solve(prob, ScipyRoot("lm")) + @test sol.retcode == ReturnCode.Success + res = system(sol.u, nothing) + @test all(abs.(res) .< 1e-10) + for method in ["broyden1", "broyden2", "anderson", "linearmixing", + "diagbroyden", "excitingmixing", "krylov", "df-sane"] + sol = solve(prob, ScipyRoot(method)) + @test sol.retcode in (ReturnCode.Success, ReturnCode.Failure) + if sol.retcode == ReturnCode.Success + res = system(sol.u, nothing) + @test all(abs.(res) .< 1e-4) + end + end + end + + @testset "ScipyLinprog" begin + function linear_obj(x, p) + c = [-1.0, -2.0] + return c + end + x0_lp = [0.0, 0.0] + optf = OptimizationFunction(linear_obj) + prob = OptimizationProblem(optf, x0_lp, nothing, + lb = [0.0, 0.0], ub = [4.0, 2.0]) + for method in ["highs", "highs-ds", "highs-ipm"] + sol = solve(prob, ScipyLinprog(method)) + @test sol.retcode in (ReturnCode.Success, ReturnCode.Failure) + if sol.retcode == ReturnCode.Success + @test sol.u[1] >= 0.0 + @test sol.u[2] >= 0.0 + @test sol.u[1] <= 4.0 + @test sol.u[2] <= 2.0 + end + end + end + + @testset "ScipyMilp" begin + function milp_obj(x, p) + c = [-1.0, -2.0] + return c + end + x0_milp = [0.0, 0.0] + optf = OptimizationFunction(milp_obj) + prob = OptimizationProblem(optf, x0_milp, nothing, + lb = [0.0, 0.0], ub = [4.0, 2.0]) + sol = solve(prob, ScipyMilp()) + @test sol.retcode in (ReturnCode.Success, ReturnCode.Failure, ReturnCode.Infeasible) + if sol.retcode == ReturnCode.Success + @test sol.u[1] >= 0.0 + @test sol.u[2] >= 0.0 + @test sol.u[1] <= 4.0 + @test sol.u[2] <= 2.0 + end + end + + @testset "cache interface" begin + objective(x, p) = (p[1] - x[1])^2 + x0 = zeros(1) + p = [1.0] + optf = OptimizationFunction(objective, Optimization.AutoZygote()) + prob = OptimizationProblem(optf, x0, p) + cache = Optimization.init(prob, ScipyBFGS()) + sol = Optimization.solve!(cache) + @test sol.retcode == ReturnCode.Success + @test sol.u ≈ [1.0] atol=1e-3 + cache = Optimization.reinit!(cache; p = [2.0]) + sol = Optimization.solve!(cache) + @test sol.u ≈ [2.0] atol=1e-3 + end + + @testset "callback" begin + cbstopping = function (state, loss) + return state.objective < 0.7 + end + optprob = OptimizationFunction(rosenbrock, Optimization.AutoZygote()) + prob = OptimizationProblem(optprob, x0, _p) + @test_throws ErrorException solve(prob, ScipyBFGS(), callback = cbstopping) + end + + @testset "constrained optimization" begin + Random.seed!(1) + cons = (res, x, p) -> res .= [x[1]^2 + x[2]^2 - 1.0] + cons_j = (res, x, p) -> begin + res[1,1] = 2*x[1] + res[1,2] = 2*x[2] + end + x0 = zeros(2) + optprob = OptimizationFunction(rosenbrock, Optimization.AutoForwardDiff(); cons = cons, cons_j = cons_j) + prob_cobyla = OptimizationProblem(optprob, x0, _p, lcons = [-1e-6], ucons = [1e-6]) + sol = solve(prob_cobyla, ScipyCOBYLA(), maxiters = 10000) + @test sol.retcode == ReturnCode.Success + @test 10 * sol.objective < l1 + Random.seed!(42) + prob = OptimizationProblem(optprob, rand(2), _p, lcons = [0.0], ucons = [0.0]) + sol = solve(prob, ScipySLSQP()) + @test sol.retcode == ReturnCode.Success + @test 10 * sol.objective < l1 + Random.seed!(123) + prob = OptimizationProblem(optprob, rand(2), _p, lcons = [0.0], ucons = [0.0]) + sol = solve(prob, ScipyTrustConstr()) + @test sol.retcode == ReturnCode.Success + @test 10 * sol.objective < l1 + function con2_c(res, x, p) + res .= [x[1]^2 + x[2]^2 - 1.0, x[2] * sin(x[1]) - x[1] - 2.0] + end + optprob = OptimizationFunction(rosenbrock, Optimization.AutoForwardDiff(); cons = con2_c) + Random.seed!(456) + prob = OptimizationProblem(optprob, rand(2), _p, lcons = [0.0, -Inf], ucons = [0.0, 0.0]) + sol = solve(prob, ScipySLSQP()) + @test sol.retcode == ReturnCode.Success + @test 10 * sol.objective < l1 + Random.seed!(789) + prob = OptimizationProblem(optprob, [0.5, 0.5], _p, lcons = [-Inf, -Inf], + ucons = [0.0, 0.0], lb = [-1.0, -1.0], ub = [1.0, 1.0]) + sol = solve(prob, ScipyShgo(), n = 50, iters = 1) + @test sol.retcode == ReturnCode.Success + @test sol.objective < l1 + end + + @testset "method-specific options" begin + simple_optprob = OptimizationFunction(rosenbrock, Optimization.AutoZygote()) + unconstrained_prob = OptimizationProblem(simple_optprob, x0, _p, lb = [-1.0, -1.0], ub = [1.0, 1.0]) + sol = solve(unconstrained_prob, ScipyDifferentialEvolution(), + popsize = 10, mutation = (0.5, 1.0), recombination = 0.7) + @test sol.retcode == ReturnCode.Success + sol = solve(unconstrained_prob, ScipyBasinhopping(), T = 1.0, stepsize = 0.5, niter = 10) + @test sol.retcode == ReturnCode.Success + sol = solve(unconstrained_prob, ScipyDualAnnealing(), + initial_temp = 5000.0, restart_temp_ratio = 2e-5) + @test sol.retcode == ReturnCode.Success + sol = solve(unconstrained_prob, ScipyShgo(), n = 50, sampling_method = "simplicial") + @test sol.retcode == ReturnCode.Success + sol = solve(unconstrained_prob, ScipyDirect(), eps = 0.001, locally_biased = true) + @test sol.retcode == ReturnCode.Success + sol = solve(unconstrained_prob, ScipyBrute(), Ns = 5, workers = 1) + @test sol.retcode == ReturnCode.Success + end + + @testset "gradient-free methods" begin + optf_no_grad = OptimizationFunction(rosenbrock) + prob = OptimizationProblem(optf_no_grad, x0, _p) + sol = solve(prob, ScipyCOBYLA(), maxiters = 10000) + @test sol.retcode == ReturnCode.Success + sol = solve(prob, ScipyNelderMead()) + @test sol.retcode == ReturnCode.Success + sol = solve(prob, ScipyPowell()) + @test sol.retcode == ReturnCode.Success + end + + @testset "AutoDiff backends" begin + for adtype in [Optimization.AutoZygote(), + Optimization.AutoReverseDiff(), + Optimization.AutoForwardDiff()] + optf = OptimizationFunction(rosenbrock, adtype) + prob = OptimizationProblem(optf, x0, _p) + sol = solve(prob, ScipyBFGS()) + @test sol.retcode == ReturnCode.Success + @test 10 * sol.objective < l1 + end + end + + @testset "optimization stats" begin + optprob = OptimizationFunction(rosenbrock, Optimization.AutoZygote()) + prob = OptimizationProblem(optprob, x0, _p) + sol = solve(prob, ScipyBFGS()) + @test sol.stats.time > 0 + end + + @testset "original result access" begin + optprob = OptimizationFunction(rosenbrock, Optimization.AutoZygote()) + prob = OptimizationProblem(optprob, x0, _p) + sol = solve(prob, ScipyBFGS()) + @test !isnothing(sol.original) + @test pyhasattr(sol.original, "success") + @test pyhasattr(sol.original, "message") + end + + @testset "tolerance settings" begin + optprob = OptimizationFunction(rosenbrock, Optimization.AutoZygote()) + prob = OptimizationProblem(optprob, x0, _p) + sol = solve(prob, ScipyNelderMead(), abstol = 1e-8) + @test sol.objective < 1e-7 + sol = solve(prob, ScipyBFGS(), reltol = 1e-8) + @test sol.objective < 1e-7 + end + + @testset "constraint satisfaction" begin + cons = (res, x, p) -> res .= [x[1]^2 + x[2]^2 - 1.0] + optprob = OptimizationFunction(rosenbrock, Optimization.AutoForwardDiff(); cons = cons) + prob = OptimizationProblem(optprob, [0.5, 0.5], _p, lcons = [-0.01], ucons = [0.01]) + sol = solve(prob, ScipySLSQP()) + @test sol.retcode == ReturnCode.Success + cons_val = [0.0] + cons(cons_val, sol.u, _p) + @test abs(cons_val[1]) < 0.011 + end + + @testset "invalid method" begin + @test_throws ArgumentError ScipyMinimize("InvalidMethodName") + @test_throws ArgumentError ScipyMinimizeScalar("InvalidMethodName") + @test_throws ArgumentError ScipyLeastSquares(method="InvalidMethodName") + @test_throws ArgumentError ScipyLeastSquares(loss="InvalidLossName") + @test_throws ArgumentError ScipyRootScalar("InvalidMethodName") + @test_throws ArgumentError ScipyRoot("InvalidMethodName") + @test_throws ArgumentError ScipyLinprog("InvalidMethodName") + end + + @testset "Edge cases" begin + f_simple(x, p) = (x[1] - p[1])^2 + prob = OptimizationProblem(f_simple, [0.0], [3.0]) + sol = solve(prob, ScipyBFGS()) + @test sol.retcode == ReturnCode.Success + @test sol.u ≈ [3.0] atol=1e-6 + optprob = OptimizationFunction(rosenbrock, Optimization.AutoZygote()) + prob = OptimizationProblem(optprob, x0, _p) + @test_throws SciMLBase.IncompatibleOptimizerError solve(prob, ScipyDifferentialEvolution()) + @test_throws SciMLBase.IncompatibleOptimizerError solve(prob, ScipyDirect()) + @test_throws SciMLBase.IncompatibleOptimizerError solve(prob, ScipyDualAnnealing()) + @test_throws SciMLBase.IncompatibleOptimizerError solve(prob, ScipyBrute()) + @test_throws ArgumentError solve(prob, ScipyBrent()) + @test_throws ArgumentError solve(prob, ScipyRootScalar("brentq")) + end + + @testset "Type stability" begin + x0_f32 = Float32[0.0, 0.0] + p_f32 = Float32[1.0, 100.0] + optprob = OptimizationFunction(rosenbrock, Optimization.AutoZygote()) + prob = OptimizationProblem(optprob, x0_f32, p_f32) + sol = solve(prob, ScipyBFGS()) + @test sol.retcode == ReturnCode.Success + @test eltype(sol.u) == Float32 + end + + @testset "ScipyLinprog matrix constraints" begin + # Minimize c^T x subject to A_ub * x <= b_ub and simple bounds + c_vec(x, p) = [1.0, 1.0] # constant cost vector + x0_lp = [0.0, 0.0] + optf_lp = OptimizationFunction(c_vec) + prob_lp = OptimizationProblem(optf_lp, x0_lp) + + A_ub = [1.0 1.0] # x1 + x2 <= 5 + b_ub = [5.0] + sol = solve(prob_lp, ScipyLinprog("highs"), + A_ub = A_ub, b_ub = b_ub, + lb = [0.0, 0.0], ub = [10.0, 10.0]) + @test sol.retcode == ReturnCode.Success + @test sol.u[1] + sol.u[2] ≤ 5.0 + 1e-8 + end + + @testset "ScipyMilp matrix constraints" begin + # Mixed-integer LP: first variable binary, second continuous + c_vec_milp(x, p) = [-1.0, -2.0] # maximize -> minimize negative + x0_milp = [0.0, 0.0] + optf_milp = OptimizationFunction(c_vec_milp) + prob_milp = OptimizationProblem(optf_milp, x0_milp) + + A = [1.0 1.0] # x1 + x2 >= 1 -> lb = 1, ub = Inf + lb_con = [1.0] + ub_con = [Inf] + integrality = [1, 0] # binary, continuous + + sol = solve(prob_milp, ScipyMilp(); + A = A, lb_con = lb_con, ub_con = ub_con, + integrality = integrality, + lb = [0.0, 0.0], ub = [1.0, 10.0]) + @test sol.retcode in (ReturnCode.Success, ReturnCode.Failure) + if sol.retcode == ReturnCode.Success + @test sol.u[1] in (0.0, 1.0) + @test isapprox(sol.u[1] + sol.u[2], 1.0; atol = 1e-6) || sol.u[1] + sol.u[2] > 1.0 + end + end +end \ No newline at end of file