diff --git a/Project.toml b/Project.toml index 135da73a..ccd4df5e 100644 --- a/Project.toml +++ b/Project.toml @@ -3,6 +3,7 @@ DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e" Ipopt = "b6b21f68-93f8-5de0-b562-5493be1d77c9" JuMP = "4076af6c-e467-56ae-b986-b466b2749572" Libdl = "8f399da3-3557-5675-b5ff-fb832c97cbdb" +LineSearches = "d3d80556-e9d4-5f37-9878-2ab0fcc64255" MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee" MsgPack = "99f44e22-a591-53d1-9472-aa23ef4bd671" NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec" diff --git a/examples/lang_julia/call_optim_rosenbrock.jl b/examples/lang_julia/call_optim_rosenbrock.jl new file mode 100644 index 00000000..bbfde2f0 --- /dev/null +++ b/examples/lang_julia/call_optim_rosenbrock.jl @@ -0,0 +1,110 @@ +using Printf + +using OpenInterfaces.Interfaces.Optim + + +function parse_args(args) + supported_impls = ["optim_jl", "ipopt_jl", "scipy_optimize"] + + impl = "optim_jl" + method = "NelderMead" + linesearch = "StrongWolfe" + + if length(args) > 0 + impl = args[1] + + if ! (impl in supported_impls) + error(""" + Given implementation $impl is not in the list\ + of the supported implementations: $supported_impls + """) + end + end + + if length(args) > 1 + method = args[2] + end + + if length(args) > 2 + linesearch = args[3] + else + linesearch + end + + return impl, method, linesearch +end + + +function rosenbrock_objective_fn(x, a) + """The Rosenbrock function with additional arguments""" + return sum(a * (x[2:end] - x[1:(end-1)] .^ 2.0) .^ 2.0 + (1 .- x[1:(end-1)]) .^ 2.0) +end + +function rosenbrock_grad_fn(x, grad_f, a) + """Gradient of the Rosenbrock function""" + + # `xi` is x[i], `xip1` is x[i+1] + xi = @view x[1:(end-1)] + xip1 = @view x[2:end] + + grad_f[1:(end-1)] .= -4.0 .* a .* xi .* (xip1 .- xi .^ 2.0) .- 2.0 .* (1.0 .- xi) + grad_f[end] = 0.0 + grad_f[2:end] .+= 2.0 .* a .* (xip1 .- xi .^ 2.0) + return 0 +end + + +function main(args) + impl, method, linesearch = parse_args(args) + println("Calling from Julia an open interface for initial-value problems (IVP)") + println("Implementation: $impl") + println("Method: $method") + println("Line search (only for optim_jl:BFGS): $linesearch") + + x0 = [3.14, 2.72, 6.18, 9.81, 8.31] + user_data = 10 + + s = Optim.Self(impl) + Optim.set_initial_guess(s, x0) + Optim.set_user_data(s, user_data) + if impl == "scipy_optimize" + if method == "NelderMead" + Optim.set_method(s, "nelder-mead", Dict("fatol" => 1e-11)) + elseif method == "BFGS" + Optim.set_method(s, "BFGS", Dict("gtol" => 1e-8)) + else + error("Unsupported method '$(method)'") + end + elseif impl == "optim_jl" + if method == "NelderMead" + Optim.set_method(s, "NelderMead", Dict("g_abstol" => 1e-11)) + elseif method == "BFGS" + Optim.set_method( + s, + "BFGS", + Dict("g_abstol" => 1e-8, "linesearch" => linesearch), + ) + else + error("Unsupported method '$(method)'") + end + else + error("Unknown implementation") + end + Optim.set_objective_fn(s, rosenbrock_objective_fn) + Optim.set_grad_fn(s, rosenbrock_grad_fn) + + status, message = Optim.minimize(s) + x = s.x + + println("Message: ", message) + @assert status == 0 + println("x = ", x) + if all(abs.(x .- 1.0) .< 1e-5) # The solution is [1, 1, ..., 1]. + println("\033[1;32mSUCCESS\033[0m Found solution is close to the exact one") + else + println("\033[1;31mFAIL\033[0m Found solution is NOT close to the exact one") + end +end + + +main(ARGS) diff --git a/examples/lang_python/call_optim_rosenbrock.py b/examples/lang_python/call_optim_rosenbrock.py new file mode 100755 index 00000000..253c2bed --- /dev/null +++ b/examples/lang_python/call_optim_rosenbrock.py @@ -0,0 +1,84 @@ +#!/usr/bin/env python3 +import argparse +import sys + +import numpy as np +from openinterfaces.interfaces.optim import Optim + + +def parse_args(): + parser = argparse.ArgumentParser() + parser.add_argument("impl", nargs="?", default="optim_jl") + parser.add_argument("method", nargs="?", default="NelderMead") + parser.add_argument("linesearch", nargs="?", default="StrongWolfe") + return parser.parse_args() + + +def rosenbrock_objective_fn(x, a): + """The Rosenbrock function with parameter""" + return np.sum(a * (x[1:] - x[:-1] ** 2.0) ** 2.0 + (1.0 - x[:-1]) ** 2.0) + + +def rosenbrock_grad_fn(x, grad_f, a): + """The gradient of the Rosenbrock function""" + xi = x[:-1] + xip1 = x[1:] + + grad_f[:-1] = -4.0 * a * xi * (xip1 - xi**2.0) - 2.0 * (1.0 - xi) + grad_f[-1] = 0.0 + grad_f[1:] += 2.0 * a * (xip1 - xi**2.0) + return 0 + + +def main(): + args = parse_args() + impl = args.impl + method = args.method + linesearch = args.linesearch + + print("Calling from Python an open interface for optimization") + print(f"Implementation: {impl}") + print(f"Method: {method}") + print(f"Line search (only for optim_jl:BFGS): {linesearch}") + + x0 = np.array([3.14, 2.72, 6.18, 9.81, 8.31]) + user_data = 10.0 + + s = Optim(impl) + s.set_initial_guess(x0) + s.set_user_data(user_data) + + if impl == "scipy_optimize": + if method == "NelderMead": + s.set_method("nelder-mead", {"fatol": 1e-11}) + elif method == "BFGS": + s.set_method("BFGS", {"gtol": 1e-8}) + else: + raise ValueError(f"Unsupported method '{method}'") + elif impl == "optim_jl": + if method == "NelderMead": + s.set_method("NelderMead", {"g_abstol": 1e-11}) + elif method == "BFGS": + s.set_method("BFGS", {"g_abstol": 1e-8, "linesearch": linesearch}) + else: + raise ValueError(f"Unsupported method '{method}'") + else: + raise ValueError("Unknown implementation") + + s.set_objective_fn(rosenbrock_objective_fn) + s.set_grad_fn(rosenbrock_grad_fn) + + status, message = s.minimize() + x = s.x + + print(f"Message: {message}") + assert status == 0 + print(f"x = {x}") + if np.all(np.abs(x - 1.0) < 1e-5): # The solution is [1, 1, ..., 1]. + print("\033[1;32mSUCCESS\033[0m Found solution is close to the exact one") + else: + print("\033[1;31mFAIL\033[0m Found solution is NOT close to the exact one") + + +if __name__ == "__main__": + sys.exit(main()) diff --git a/lang_julia/OpenInterfacesImpl/src/serialization.jl b/lang_julia/OpenInterfacesImpl/src/serialization.jl index cac25523..c8d727e6 100644 --- a/lang_julia/OpenInterfacesImpl/src/serialization.jl +++ b/lang_julia/OpenInterfacesImpl/src/serialization.jl @@ -22,6 +22,8 @@ function deserialize(sd::Ptr, len::Csize_t)::Dict # deliberately to wider integer type. if typeof(elem) == UInt8 elem = Int64(elem) + elseif typeof(elem) == UInt16 + elem = Int64(elem) end if i % 2 == 1 push!(data, Symbol(elem)) diff --git a/lang_julia/_impl/optim/optim_jl/optim_jl.jl b/lang_julia/_impl/optim/optim_jl/optim_jl.jl index 5df080a0..153676a3 100644 --- a/lang_julia/_impl/optim/optim_jl/optim_jl.jl +++ b/lang_julia/_impl/optim/optim_jl/optim_jl.jl @@ -4,14 +4,18 @@ module OptimJl export Self, set_initial_guess +using LineSearches using Optim GENERAL_OPTIONS_NAMES = [ :x_abstol, :x_reltol, + :x_tol, # Deprecated, keep to be user-friendly :f_abstol, :f_reltol, + :f_tol, # Deprecated, keep to be user-friendly :g_abstol, + :g_tol, # Deprecated, keep to be user-friendly :f_calls_limit, :g_calls_limit, :h_calls_limit, @@ -29,12 +33,11 @@ mutable struct Self grad_fn::Union{Function,Nothing} user_data::Any method_name::String - method_params::Dict + method_options::Dict method::Any general_options::Dict - grad_values::Any function Self() - return new([], nothing, nothing, nothing, "BFGS", Dict(), BFGS(), Dict(), nothing) + return new([], nothing, nothing, nothing, "BFGS", Dict(), BFGS(), Dict()) end end @@ -54,7 +57,7 @@ end function set_grad_fn(self::Self, grad_fn) self.grad_fn = grad_fn - self.grad_values = zeros(length(self.x0)) + self.grad_fn(self.x0, zeros(length(self.x0)), self.user_data) end function set_method(self, method_name, method_params) @@ -65,6 +68,12 @@ function set_method(self, method_name, method_params) general_options = Dict() method_options = Dict() + if haskey(method_params, :linesearch) + method_options[:linesearch] = + getfield(LineSearches, Symbol(method_params[:linesearch]))() + delete!(method_params, :linesearch) + end + for (k, v) in method_params if k in GENERAL_OPTIONS_NAMES general_options[k] = v @@ -78,9 +87,11 @@ function set_method(self, method_name, method_params) error("Unknown or unsupported method '$method_name'") end self.method_name = method_name - self.method_params = merge(self.method_params, method_params) + self.method_options = merge(self.method_options, method_options) self.method = getfield(Optim, method_symbol)(; method_options...) self.general_options = general_options + println("method_options = ", method_options) + println("general_options = ", general_options) end function minimize(self::Self, out_x)::Tuple{Int,String} diff --git a/lang_python/oif_impl/openinterfaces/_impl/optim/scipy_optimize/scipy_optimize.py b/lang_python/oif_impl/openinterfaces/_impl/optim/scipy_optimize/scipy_optimize.py index 13163edd..f7c971ad 100644 --- a/lang_python/oif_impl/openinterfaces/_impl/optim/scipy_optimize/scipy_optimize.py +++ b/lang_python/oif_impl/openinterfaces/_impl/optim/scipy_optimize/scipy_optimize.py @@ -94,6 +94,8 @@ def minimize(self, out_x): out_x[:] = result.x print(result) + if hasattr(result, "jac"): + print("|∇f(x*)|_∞ = ", np.linalg.norm(result.jac, np.inf)) return (result.status, result.message) def grad_wrapper(self, x, __): diff --git a/tests/lang_julia/test_examples.jl b/tests/lang_julia/test_examples.jl index caffd4bc..f4ea5da7 100644 --- a/tests/lang_julia/test_examples.jl +++ b/tests/lang_julia/test_examples.jl @@ -6,6 +6,7 @@ const CALL_LINSOLVE = joinpath(EXAMPLES_PATH, "call_linsolve.jl") const CALL_IVP = joinpath(EXAMPLES_PATH, "call_ivp.jl") const CALL_IVP_BURGERS = joinpath(EXAMPLES_PATH, "call_ivp_burgers_eq.jl") const CALL_IVP_VDP = joinpath(EXAMPLES_PATH, "call_ivp_vdp.jl") +const CALL_OPTIM_ROSENBROCK = joinpath(EXAMPLES_PATH, "call_optim_rosenbrock.jl") @testset "Testing Julia examples" begin @testset "Testing QEQ example with all implementations" begin @@ -138,4 +139,28 @@ const CALL_IVP_VDP = joinpath(EXAMPLES_PATH, "call_ivp_vdp.jl") end # END Van der Pol equation: unsuccessful and successful runs. # -------------------------------------------------------------------------- + + # -------------------------------------------------------------------------- + # BEGIN Optimization example with the Rosenbrock function. + @testset "Testing Optim Rosenbrock example with all implementations" begin + for impl in ["scipy_optimize", "optim_jl"] + for method in ["NelderMead", "BFGS"] + @testset "call_optim_rosenbrock.jl $impl fails due to numerics" begin + io_err = IOBuffer() + process = try + run(pipeline(`julia $CALL_OPTIM_ROSENBROCK $impl`, stderr = io_err)) + catch e + end + captured_stderr = String(take!(io_err)) + println("Captured stderr:") + println("--- BEGIN") + println(captured_stderr) + println("--- END") + @test process.exitcode == 0 + end + end + end + end + # END Optimization example with the Rosenbrock function. + # -------------------------------------------------------------------------- end diff --git a/tests/lang_python/test_examples.py b/tests/lang_python/test_examples.py index 84f062c5..f963f7c1 100644 --- a/tests/lang_python/test_examples.py +++ b/tests/lang_python/test_examples.py @@ -92,3 +92,30 @@ def test_call_ivp_from_python_vdp__with_successful_impl_succeeds(successful_impl text=True, ) assert p.returncode == 0 + + +# ----------------------------------------------------------------------------- +# Examples of optimization with the Rosenbrock function. +@pytest.fixture(params=["NelderMead", "BFGS"]) +def optim_method(request): + return request.param + + +@pytest.fixture(params=["scipy_optimize", "optim_jl"]) +def optim_impl(request): + return request.param + + +def test_call_optim_rosenbrock_from_python__exit_success( + optim_impl: str, optim_method: str +): + p = subprocess.run( + [ + "python", + "examples/lang_python/call_optim_rosenbrock.py", + optim_impl, + optim_method, + ] + ) + + assert p.returncode == 0