diff --git a/.github/workflows/qa.yaml b/.github/workflows/qa.yaml index 40f8eefc7..a4939fe0e 100644 --- a/.github/workflows/qa.yaml +++ b/.github/workflows/qa.yaml @@ -14,7 +14,7 @@ jobs: strategy: matrix: python-version: ['3.10', '3.11', '3.12', '3.13', '3.14'] - julia-version: ['1.11', '1.12'] + julia-version: ['1.11'] # Excluding Julia 1.12 due to memory crashes , '1.12'] defaults: run: shell: bash -l {0} diff --git a/Makefile b/Makefile index 1c81c3dab..ec4bb4c4e 100644 --- a/Makefile +++ b/Makefile @@ -32,7 +32,7 @@ help: .PHONY : test test : @echo "=== C tests ===" - cd build && ctest --output-on-failure + cd build && ctest --verbose --output-on-failure @echo "=== Julia tests ===" julia tests/lang_julia/runtests.jl @echo "=== Python tests ===" diff --git a/Project.toml b/Project.toml index 44d7156d2..135da73a9 100644 --- a/Project.toml +++ b/Project.toml @@ -1,4 +1,5 @@ [deps] +DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e" Ipopt = "b6b21f68-93f8-5de0-b562-5493be1d77c9" JuMP = "4076af6c-e467-56ae-b986-b466b2749572" Libdl = "8f399da3-3557-5675-b5ff-fb832c97cbdb" @@ -16,3 +17,9 @@ SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462" [sources] OpenInterfaces = {path = "lang_julia/OpenInterfaces"} OpenInterfacesImpl = {path = "lang_julia/OpenInterfacesImpl"} + +[compat] +MsgPack = "1.2.1" +Optim = "2.0.1" +OrdinaryDiffEq = "6.102.1" +SciMLBase = "2.120.0" diff --git a/lang_julia/OpenInterfaces/src/OpenInterfaces.jl b/lang_julia/OpenInterfaces/src/OpenInterfaces.jl index 31194938b..097152b36 100644 --- a/lang_julia/OpenInterfaces/src/OpenInterfaces.jl +++ b/lang_julia/OpenInterfaces/src/OpenInterfaces.jl @@ -10,6 +10,7 @@ export OIF_TYPE_INT, OIF_TYPE_STRING, OIF_TYPE_CALLBACK, OIF_USER_DATA, + OIF_TYPE_USER_DATA, OIF_TYPE_CONFIG_DICT # Language ids @@ -100,6 +101,7 @@ mutable struct OIFArrayF64 flags::Int32 function OIFArrayF64(arr::AbstractArray{Float64}) + nd = ndims(arr) dims = collect(Cssize_t, size(arr)) flags = OIF_ARRAY_F_CONTIGUOUS @@ -107,14 +109,18 @@ mutable struct OIFArrayF64 flags = OIF_ARRAY_C_CONTIGUOUS | OIF_ARRAY_F_CONTIGUOUS end - self = new( - ndims(arr), - Base.unsafe_convert(Ptr{Cssize_t}, dims), - Base.unsafe_convert(Ptr{Float64}, arr), - flags, - ) + dimensions = Ptr{Cssize_t}(Libc.malloc(nd * sizeof(Cssize_t))) + unsafe_copyto!(dimensions, pointer(dims), nd) + + self = new(nd, dimensions, Base.unsafe_convert(Ptr{Float64}, arr), flags) + + finalizer(finalizing, self) - return GC.@preserve dims self + return self + end + + function finalizing(self) + Libc.free(self.dimensions) end end @@ -558,7 +564,7 @@ end # otherwise, the will be crashes or something like that # as Julia's garbage collector is pretty fast to remove reference # that it thinks are unused. -function make_oif_user_data(data::Ref{Any})::OIFUserData +function make_oif_user_data(data)::OIFUserData data_ptr = Base.unsafe_convert(Ptr{Cvoid}, data) OIFUserData(OIF_LANG_JULIA, C_NULL, data_ptr, C_NULL) end diff --git a/lang_julia/OpenInterfacesImpl/Project.toml b/lang_julia/OpenInterfacesImpl/Project.toml index 7a0191ff5..929fc6c85 100644 --- a/lang_julia/OpenInterfacesImpl/Project.toml +++ b/lang_julia/OpenInterfacesImpl/Project.toml @@ -1,7 +1,7 @@ name = "OpenInterfacesImpl" uuid = "daff4927-7858-4806-96da-915047fb441c" -authors = ["Dmitry Kabanov ", "Stephan Rave ", "Stephan Rave begin + if argtype == OIF_TYPE_I32 + Cint + elseif argtype == OIF_TYPE_F64 + Cdouble + elseif argtype == OIF_TYPE_ARRAY_F64 + Ref{OIFArrayF64} + elseif argtype == OIF_TYPE_USER_DATA + Ptr{Cvoid} + else + error("Unsupported argument type: $argtype") + end + end, oif_argtypes), + ) -function make_wrapper_over_c_callback(fn_c::Ptr{Cvoid})::Function - function wrapper(t::Float64, y::Vector{Float64}, ydot::Vector{Float64}, user_data)::Int - if typeof(user_data) == SciMLBase.NullParameters - user_data = C_NULL - end + # Build a concrete Tuple type e.g. Tuple{Cdouble,Cdouble}. + # This IS a valid type parameter, unlike a plain tuple value (Cdouble, Cdouble). + c_argtypes_type = Core.apply_type(Tuple, c_argtypes...) - oif_y = Ref(OIFArrayF64(y)) - oif_ydot = Ref(OIFArrayF64(ydot)) + c_restype = if oif_restype == OIF_TYPE_I32 || oif_restype == OIF_TYPE_INT + # The second condition is just for explicitness, as `INT` is an alias to `I32`. + Cint + elseif oif_restype == OIF_TYPE_F64 + Cdouble + else + error( + "Unsupported return type: $oif_restype. Only supported types are Cint and Cdouble", + ) + end - @ccall $fn_c( - t::Float64, - oif_y::Ptr{OIFArrayF64}, - oif_ydot::Ptr{OIFArrayF64}, - user_data::Ptr{Cvoid}, - )::Cint + function wrapper(args...) + c_args_from_jl_args = Tuple( + map( + (argtype, argvalue) -> begin + if argtype == OIF_TYPE_I32 + argvalue + elseif argtype == OIF_TYPE_F64 + argvalue + elseif argtype == OIF_TYPE_ARRAY_F64 + OIFArrayF64(argvalue) + elseif argtype == OIF_TYPE_USER_DATA + argvalue === nothing ? C_NULL : argvalue + else + error("Unsupported argument type: $argtype") + end + end, + oif_argtypes, + args, + ), + ) - return 0 + return GC.@preserve args c_args_from_jl_args call_c( + fn_c, + c_argtypes_type, + c_restype, + c_args_from_jl_args, + ) end + return wrapper end +@generated function call_c( + fn_p, + ::Type{argtypes}, + ::Type{restype}, + args, +) where {argtypes,restype} + + n = length(argtypes.parameters) + + argtypes_expr = Expr(:tuple, map(t -> :($t), argtypes.parameters)...) + # Generate :(args[1]), :(args[2]), … at code-generation time so that ccall + # sees a fixed, statically-known argument list (runtime splatting is not allowed). + args_expr = [:(args[$i]) for i = 1:n] + + quote + ccall(fn_p, $restype, $argtypes_expr, $(args_expr...)) + end +end + end diff --git a/lang_julia/_impl/ivp/jl_diffeq/jl_diffeq.jl b/lang_julia/_impl/ivp/jl_diffeq/jl_diffeq.jl index df2edfd5f..caeb57a8b 100644 --- a/lang_julia/_impl/ivp/jl_diffeq/jl_diffeq.jl +++ b/lang_julia/_impl/ivp/jl_diffeq/jl_diffeq.jl @@ -3,6 +3,7 @@ export Self, set_initial_value, set_rhs_fn, set_tolerances, integrate, set_user_data, set_integrator using OrdinaryDiffEq +import SciMLBase mutable struct Self t0::Float64 @@ -79,6 +80,10 @@ function _rhs_wrapper(rhs) # Note that the `p` argument will always be in the function signature # regardless of if the problem is defined with parameters! function wrapper(du, u, p, t) + if typeof(p) == SciMLBase.NullParameters + return rhs(t, u, du, nothing) + end + return rhs(t, u, du, p) end return wrapper diff --git a/lang_julia/_impl/optim/optim_jl/optim_jl.jl b/lang_julia/_impl/optim/optim_jl/optim_jl.jl index a27128ecd..efbaca5e2 100644 --- a/lang_julia/_impl/optim/optim_jl/optim_jl.jl +++ b/lang_julia/_impl/optim/optim_jl/optim_jl.jl @@ -47,6 +47,7 @@ end function set_objective_fn(self::Self, objective_fn) self.objective_fn = objective_fn + self.objective_fn(self.x0, self.user_data) end function set_method(self, method_name, method_params) diff --git a/lang_julia/bridge_julia.c b/lang_julia/bridge_julia.c index cc87af349..e4d6b8b05 100644 --- a/lang_julia/bridge_julia.c +++ b/lang_julia/bridge_julia.c @@ -203,10 +203,85 @@ build_julia_tuple_from_size_t_array(intptr_t *dimensions, size_t ndims) return tuple; } +/** + * Build a Julia tuple from an `int32_t *` array. + * + * @param vals Array with elements of type int (32-bit integers) + * @param n Number of the elements in the array + * @return Julia tuple on success, `NULL` otherwise + */ +jl_value_t * +build_julia_tuple_from_oifargtype_array(OIFArgType *vals, size_t n) +{ + char *tuple_string = malloc(sizeof(char) * 1024); + char *cursor = tuple_string; + char *right_bound = &tuple_string[1024 - 1]; + intptr_t diff = right_bound - cursor; + snprintf(tuple_string, diff, "%s", "("); + cursor++; + char buffer[BUFFER_SIZE_]; + int chars_written; + jl_value_t *tuple = NULL; + + for (size_t i = 0; i < n; ++i) { + chars_written = snprintf(buffer, BUFFER_SIZE_, "%d", vals[i]); + assert(chars_written >= 1); + if (chars_written >= BUFFER_SIZE_) { + goto report_too_long; + } + diff = right_bound - cursor; + chars_written = snprintf(cursor, diff, "%s", buffer); + assert(chars_written >= 1); + if (chars_written >= diff) { + goto report_too_long; + } + cursor += chars_written; + diff = right_bound - cursor; + if (i < n - 1 || n == 1) { + chars_written = snprintf(cursor, diff, "%s", ","); + assert(chars_written >= 1); + if (chars_written >= diff) { + goto report_too_long; + } + cursor++; + } + } + diff = right_bound - cursor; + chars_written = snprintf(cursor, diff, "%s", ")"); + if (chars_written >= diff) { + fprintf(stderr, "ERROR: the string to copy is too long\n"); + exit(1); + } + + tuple = jl_eval_string(tuple_string); + if (jl_exception_occurred()) { + const char *p = jl_string_ptr( + jl_eval_string("sprint(showerror, ccall(:jl_exception_occurred, Any, ()))")); + fprintf(stderr, "%s\n", p); + } + goto cleanup; + +report_too_long: + fprintf(stderr, + "[build_julia_tuple_from_int32_t_array] The string representation of the julia " + "tuple does not fit in the buffer\n"); + +cleanup: + free(tuple_string); + + return tuple; +} + static jl_value_t * make_wrapper_over_c_callback(OIFCallback *p) { jl_value_t *wrapper = NULL; + + jl_value_t *fn_p_c_wrapped = NULL; + jl_value_t *oif_argtypes = NULL; + jl_value_t *oif_restype = NULL; + JL_GC_PUSH3(&fn_p_c_wrapped, &oif_argtypes, &oif_restype); // NOLINT + if (CALLBACK_MODULE_ == NULL) { jl_eval_string("using OpenInterfacesImpl.CallbackWrapper"); if (jl_exception_occurred()) { @@ -224,10 +299,14 @@ make_wrapper_over_c_callback(OIFCallback *p) jl_get_function(CALLBACK_MODULE_, "make_wrapper_over_c_callback"); assert(fn_callback != NULL); assert(p->fn_p_c != NULL); - jl_value_t *fn_p_c_wrapped = jl_box_voidpointer(p->fn_p_c); - wrapper = jl_call1(fn_callback, fn_p_c_wrapped); + fn_p_c_wrapped = jl_box_voidpointer(p->fn_p_c); + oif_argtypes = build_julia_tuple_from_oifargtype_array(p->arg_types, p->nargs); + oif_restype = jl_box_int32(p->restype); + + wrapper = jl_call3(fn_callback, fn_p_c_wrapped, oif_argtypes, oif_restype); cleanup: + JL_GC_POP(); return wrapper; } diff --git a/lang_python/oif/openinterfaces/core.py b/lang_python/oif/openinterfaces/core.py index 657ea444d..6e37e5f34 100644 --- a/lang_python/oif/openinterfaces/core.py +++ b/lang_python/oif/openinterfaces/core.py @@ -139,6 +139,10 @@ class OIFCallback(ctypes.Structure): ("fn_p_c", ctypes.c_void_p), ("fn_p_jl", ctypes.c_void_p), ("fn_p_py", ctypes.c_void_p), + ("nargs", ctypes.c_uint32), + ("arg_types", ctypes.POINTER(OIFArgType)), + ("restype", OIFArgType), + ("hidden", ctypes.c_void_p), ] @@ -197,7 +201,15 @@ def make_oif_callback( fn_p_c = ctypes.cast(c_wrapper_fn, ctypes.c_void_p) assert fn_p_c.value is not None - oifcallback = OIFCallback(OIF_LANG_PYTHON, fn_p_c, None, fn_p_py) + + nargs = len(argtypes) + oif_arg_types = ctypes.cast( + (ctypes.c_int * nargs)(*argtypes), ctypes.POINTER(OIFArgType) + ) + + oifcallback = OIFCallback( + OIF_LANG_PYTHON, fn_p_c, None, fn_p_py, nargs, oif_arg_types, restype, None + ) return oifcallback 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 fe8b4a9f6..b4d939e0c 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 @@ -84,6 +84,7 @@ def minimize(self, out_x): ) out_x[:] = result.x + print(result) return (result.status, result.message) # def print_stats(self): diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 4a99049e4..8111ae41a 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -25,3 +25,4 @@ endif() include(GoogleTest) add_subdirectory(lang_c) +add_subdirectory(lang_julia) diff --git a/tests/lang_julia/CMakeLists.txt b/tests/lang_julia/CMakeLists.txt new file mode 100644 index 000000000..2b5ee7cca --- /dev/null +++ b/tests/lang_julia/CMakeLists.txt @@ -0,0 +1,2 @@ +add_library(oif_aux_funcs SHARED aux_funcs.c) +target_include_directories(oif_aux_funcs PRIVATE ${CMAKE_SOURCE_DIR}/include) diff --git a/tests/lang_julia/aux_funcs.c b/tests/lang_julia/aux_funcs.c new file mode 100644 index 000000000..b2de78f65 --- /dev/null +++ b/tests/lang_julia/aux_funcs.c @@ -0,0 +1,54 @@ +#include + +typedef struct context { + int p1; + double p2; +} Context; + +int +add_two_ints(int a, int b) +{ + return a + b; +} + +double +add_two_doubles(double a, double b) +{ + return a + b; +} + +int +add_i32_f64__ret_i32(int a, double b) +{ + return (int)(a + b); +} + +double +add_i32_f64__ret_f64(int a, double b) +{ + return a + b; +} + +int +axpy_f64_arrf64_arrf64__ret_i32(double alpha, OIFArrayF64 *x, OIFArrayF64 *y) +{ + int N = x->dimensions[0]; + for (int i = 0; i < N; i++) { + y->data[i] = y->data[i] + alpha * x->data[i]; + } + + return 0; +} + +int +axpy_with_context_f64_arrf64_arrf64__ret_i32(double alpha, OIFArrayF64 *x, OIFArrayF64 *y, + void *void_context) +{ + Context *context = void_context; + int N = x->dimensions[0]; + for (int i = 0; i < N; i++) { + y->data[i] = y->data[i] + (alpha + context->p2) * (x->data[i] + context->p1); + } + + return 0; +} diff --git a/tests/lang_julia/test_callback.jl b/tests/lang_julia/test_callback.jl new file mode 100644 index 000000000..3eb5fd5f2 --- /dev/null +++ b/tests/lang_julia/test_callback.jl @@ -0,0 +1,74 @@ +using Test + +using OpenInterfaces +using OpenInterfacesImpl.CallbackWrapper +using Libdl + + +mutable struct Context + p1::Int32 + p2::Float64 +end + +@testset "check that callback wrapper can be dynamic" begin + lib = Libdl.dlopen("liboif_aux_funcs.so") + add_two_ints_fn = Libdl.dlsym(lib, :add_two_ints) + add_two_doubles_fn = Libdl.dlsym(lib, :add_two_doubles) + add_i32_f64__ret_i32_fn = Libdl.dlsym(lib, :add_i32_f64__ret_i32) + add_i32_f64__ret_f64_fn = Libdl.dlsym(lib, :add_i32_f64__ret_f64) + axpy_f64_arrf64_arrf64__ret_i32_fn = Libdl.dlsym(lib, :axpy_f64_arrf64_arrf64__ret_i32) + axpy_with_context_f64_arrf64_arrf64__ret_i32_fn = + Libdl.dlsym(lib, :axpy_with_context_f64_arrf64_arrf64__ret_i32) + + w_f64_f64__ret_f64 = make_wrapper_over_c_callback( + add_two_doubles_fn, + (OIF_TYPE_F64, OIF_TYPE_F64), + OIF_TYPE_F64, + ) + @test w_f64_f64__ret_f64(1.4, 2.7) ≈ 4.1 + + w_i32_i32__ret_i32 = make_wrapper_over_c_callback( + add_two_ints_fn, + (OIF_TYPE_I32, OIF_TYPE_I32), + OIF_TYPE_I32, + ) + @test w_i32_i32__ret_i32(42, 69) == 111 + + add_i32_f64__ret_i32_wrapper = make_wrapper_over_c_callback( + add_i32_f64__ret_i32_fn, + (OIF_TYPE_I32, OIF_TYPE_F64), + OIF_TYPE_I32, + ) + @test add_i32_f64__ret_i32_wrapper(12, 3.14) == 15 + + add_i32_f64__ret_f64_wrapper = make_wrapper_over_c_callback( + add_i32_f64__ret_i32_fn, + (OIF_TYPE_I32, OIF_TYPE_F64), + OIF_TYPE_F64, + ) + @test add_i32_f64__ret_f64_wrapper(12, 3.14) ≈ 15.14 + + axpy_f64_arrf64_arrf64__ret_i32_wrapper = make_wrapper_over_c_callback( + axpy_f64_arrf64_arrf64__ret_i32_fn, + (OIF_TYPE_F64, OIF_TYPE_ARRAY_F64, OIF_TYPE_ARRAY_F64), + OIF_TYPE_I32, + ) + alpha = 2.3 + x = [1.2, 3.4, 5.6] + y = [27.25, 91.50, 16.16] + axpy_f64_arrf64_arrf64__ret_i32_wrapper(alpha, x, y) + @test y ≈ [30.01, 99.32, 29.04] + + axpy_with_context_f64_arrf64_arrf64__ret_i32_wrapper = make_wrapper_over_c_callback( + axpy_with_context_f64_arrf64_arrf64__ret_i32_fn, + (OIF_TYPE_F64, OIF_TYPE_ARRAY_F64, OIF_TYPE_ARRAY_F64, OIF_TYPE_USER_DATA), + OIF_TYPE_I32, + ) + alpha = 2.3 + x = [1.2, 3.4, 5.6] + y = [27.25, 91.50, 16.16] + context = Context(42, 3.14) + context_ref = Ref(context) + axpy_with_context_f64_arrf64_arrf64__ret_i32_wrapper(alpha, x, y, context_ref) + @test y ≈ [262.258, 338.476, 275.104] +end diff --git a/tests/lang_python/test_optim.py b/tests/lang_python/test_optim.py index 5cbb09beb..8c79442fb 100644 --- a/tests/lang_python/test_optim.py +++ b/tests/lang_python/test_optim.py @@ -2,6 +2,21 @@ import pytest from openinterfaces.interfaces.optim import Optim +# ----------------------------------------------------------------------------- +# Finding implementations. +# POSSIBLE_IMPLEMENTATIONS = ["scipy_optimize", "ipopt_jl", "optim_jl"] +POSSIBLE_IMPLEMENTATIONS = ["scipy_optimize", "optim_jl"] +IMPLEMENTATIONS = [] + +for impl in POSSIBLE_IMPLEMENTATIONS: + try: + Optim(impl) + IMPLEMENTATIONS.append(impl) + except ValueError as e: + print(f"WARNING: {str(e)}") + print(f"WARNING: Skipping tests with implementation {impl}") +# ----------------------------------------------------------------------------- + # ----------------------------------------------------------------------------- # Problems @@ -19,13 +34,13 @@ def rosenbrock_objective_fn(x, __): return sum(a * (x[1:] - x[:-1] ** 2.0) ** 2.0 + (1 - x[:-1]) ** 2.0) + b -@pytest.fixture(params=["scipy_optimize"]) +@pytest.fixture(params=IMPLEMENTATIONS) def s(request): return Optim(request.param) # ----------------------------------------------------------------------------- -# Tests +# General tests def test__simple_convex_problem__converges(s): x0 = np.array([0.5, 0.6, 0.7]) @@ -71,25 +86,32 @@ def test__parameterized_convex_problem__converges(s): assert np.all(np.abs(x - user_data) < 1e-6) -def test__parameterized_convex_problem__converges_better_with_tigher_tolerance(s): - x0 = np.array([0.5, 0.6, 0.7]) - user_data = np.array([2.0, 7.0, -1.0]) +def test__set_method_with__wrong_method_params__are_not_accepted(s): + with pytest.raises(RuntimeError): + s.set_method("BFGS", {"wrong_param": 42}) - s.set_initial_guess(x0) - s.set_user_data(user_data) - s.set_objective_fn(convex_objective_with_args_fn) - s.set_method("nelder-mead", {"xatol": 1e-6}) - status, message = s.minimize() - x_1 = s.x.copy() +# END General tests +# ----------------------------------------------------------------------------- - s.set_method("nelder-mead", {"xatol": 1e-8}) - status, message = s.minimize() - x_2 = s.x.copy() - assert np.linalg.norm(x_2 - user_data) < np.linalg.norm(x_1 - user_data) +if "scipy_optimize" in IMPLEMENTATIONS: + def test__parameterized_convex_problem__converges_better_with_tigher_tolerance(): + x0 = np.array([0.5, 0.6, 0.7]) + user_data = np.array([2.0, 7.0, -1.0]) -def test__set_method_with__wrong_method_params__are_not_accepted(s): - with pytest.raises(RuntimeError): - s.set_method("nelder-mead", {"wrong_param": 42}) + s = Optim("scipy_optimize") + s.set_initial_guess(x0) + s.set_user_data(user_data) + s.set_objective_fn(convex_objective_with_args_fn) + + s.set_method("nelder-mead", {"xatol": 1e-6}) + status, message = s.minimize() + x_1 = s.x.copy() + + s.set_method("nelder-mead", {"xatol": 1e-8}) + status, message = s.minimize() + x_2 = s.x.copy() + + assert np.linalg.norm(x_2 - user_data) < np.linalg.norm(x_1 - user_data)