Updates to the optim interface and supported implementations

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"

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)

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())

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))

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}

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, __):

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

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