diff --git a/.clang-tidy b/.clang-tidy index 9fcad0e59..61621e54c 100644 --- a/.clang-tidy +++ b/.clang-tidy @@ -1,5 +1,4 @@ --- -# This file is taken from: https://nrk.neocities.org/articles/c-static-analyzers Checks: > performance-*, misc-*, @@ -26,3 +25,5 @@ WarningsAsErrors: '*' CheckOptions: - key: bugprone-assert-side-effect.AssertMacros value: 'ASSERT' + - key: misc-include-cleaner.IgnoreHeaders + value: '[""]' diff --git a/dispatch/dispatch.c b/dispatch/dispatch.c index 404410360..eac89e88e 100644 --- a/dispatch/dispatch.c +++ b/dispatch/dispatch.c @@ -321,7 +321,8 @@ unload_interface_impl(ImplHandle implh) } int -call_interface_impl(ImplHandle implh, const char *method, OIFArgs *in_args, OIFArgs *out_args) +call_interface_impl(ImplHandle implh, const char *method, OIFArgs *in_args, OIFArgs *out_args, + OIFArgs *return_args) { int status; @@ -343,7 +344,7 @@ call_interface_impl(ImplHandle implh, const char *method, OIFArgs *in_args, OIFA } void *lib_handle = OIF_DISPATCH_HANDLES[dh]; - int (*call_impl_fn)(ImplInfo *, const char *, OIFArgs *, OIFArgs *); + int (*call_impl_fn)(ImplInfo *, const char *, OIFArgs *, OIFArgs *, OIFArgs *); call_impl_fn = dlsym(lib_handle, "call_impl"); if (call_impl_fn == NULL) { logerr(prefix_, @@ -352,7 +353,7 @@ call_interface_impl(ImplHandle implh, const char *method, OIFArgs *in_args, OIFA OIF_LANG_FROM_LANG_ID[dh]); return -3; } - status = call_impl_fn(impl_info, method, in_args, out_args); + status = call_impl_fn(impl_info, method, in_args, out_args, return_args); if (status) { logerr(prefix_, diff --git a/docs/source/api/api-python/optim/index.rst b/docs/source/api/api-python/optim/index.rst new file mode 100644 index 000000000..7e847d53e --- /dev/null +++ b/docs/source/api/api-python/optim/index.rst @@ -0,0 +1,135 @@ +optim +===== + +.. py:module:: optim + +.. autoapi-nested-parse:: + + This module defines the interface for solving minimization problems: + + .. math:: + minimize_x f(x) + + where :math:`f : \mathbb R^n \to \mathbb R`. + + + + + + + +Module Contents +--------------- + +.. py:type:: ObjectiveFn + :canonical: Callable[[np.ndarray, object], int] + + + Signature of the objective function :math:`f(, y)`. + + !!!!!!!!! The function accepts four arguments: + !!!!!!!!! - `t`: current time, + !!!!!!!!! - `y`: state vector at time :math:`t`, + !!!!!!!!! - `ydot`: output array to which the result of function evalutation is stored, + !!!!!!!!! - `user_data`: additional context (user-defined data) that + !!!!!!!!! must be passed to the function (e.g., parameters of the system). + +.. py:class:: OptimResult + + .. py:attribute:: status + :type: int + + + .. py:attribute:: x + :type: numpy.ndarray + + +.. py:class:: Optim(impl: str) + + Interface for solving optimization (minimization) problems. + + This class serves as a gateway to the implementations of the + solvers for optimization problems. + + :param impl: Name of the desired implementation. + :type impl: str + + .. rubric:: Examples + + Let's solve the following initial value problem: + + .. math:: + y'(t) = -y(t), \quad y(0) = 1. + + First, import the necessary modules: + >>> import numpy as np + >>> from oif.interfaces.ivp import IVP + + Define the right-hand side function: + + >>> def rhs(t, y, ydot, user_data): + ... ydot[0] = -y[0] + ... return 0 # No errors, optional + + Now define the initial condition: + + >>> y0, t0 = np.array([1.0]), 0.0 + + Create an instance of the IVP solver using the implementation "jl_diffeq", + which is an adapter to the `OrdinaryDiffeq.jl` Julia package: + + >>> s = IVP("jl_diffeq") + + We set the initial value, the right-hand side function, and the tolerance: + + >>> s.set_initial_value(y0, t0) + >>> s.set_rhs_fn(rhs) + >>> s.set_tolerances(1e-6, 1e-12) + + Now we integrate to time `t = 1.0` in a loop, outputting the current value + of `y` with time step `0.1`: + + >>> t = t0 + >>> times = np.linspace(t0, t0 + 1.0, num=11) + >>> for t in times[1:]: + ... s.integrate(t) + ... print(f"{t:.1f} {s.y[0]:.6f}") + 0.1 0.904837 + 0.2 0.818731 + 0.3 0.740818 + 0.4 0.670320 + 0.5 0.606531 + 0.6 0.548812 + 0.7 0.496585 + 0.8 0.449329 + 0.9 0.406570 + 1.0 0.367879 + + + .. py:attribute:: x0 + :type: numpy.ndarray + + Current value of the state vector. + + + .. py:attribute:: status + :value: -1 + + + + .. py:attribute:: x + :type: numpy.ndarray + + + .. py:method:: set_initial_guess(x0: numpy.ndarray) + + Set initial guess for the optimization problem + + + + .. py:method:: set_objective_fn(objective_fn: ObjectiveFn) + + + .. py:method:: minimize() + + Integrate to time `t` and write solution to `y`. diff --git a/docs/source/technotes/2025-05-05-how-to-release.md b/docs/source/technotes/2025-05-05-how-to-release.md index 27b4ed90b..4125db7b3 100644 --- a/docs/source/technotes/2025-05-05-how-to-release.md +++ b/docs/source/technotes/2025-05-05-how-to-release.md @@ -5,7 +5,10 @@ * `pyproject.toml` * `version.py` * `conf.py` -- [ ] Run pre-commit + * `Project.toml` + * `Manifest.toml` + * `CITATION.cff` +- [ ] Run `pre-commit run --all-files` - [ ] Make a pull request and make sure it passes all the checks - [ ] Merge the pull request to `main` - [ ] Go to the GitHub page of the repository and press "Releases" diff --git a/docs/source/technotes/2025-11-20-optimization-interface.md b/docs/source/technotes/2025-11-20-optimization-interface.md new file mode 100644 index 000000000..523123243 --- /dev/null +++ b/docs/source/technotes/2025-11-20-optimization-interface.md @@ -0,0 +1,150 @@ +# 2025-11-20 Optimization interface + +Optimization interface should accommodate adapting implementations +for general constrained nonlinear programming: +```math + minimize_{x \in \mathbb R^n} &f(x) + subject to & l <= x <= u, + & l2 <= g(x) <= u2, +``` +where $f(x) : \mathbb R^n \to \mathbb R$ and $g(x)$ are nonlinear functions of $x$. + +## Julia package JuMP.jl + +Julia package /JuMP/ (Julia Mathematical Programming) is a metapackage +that contains a lot of packages for different types of optimization problems: +for example, for integer programming, for linear programming, nonlinear, etc. +It also supports solvers with commercial licenses. + +See: +https://jump.dev/JuMP.jl/stable/tutorials/getting_started/getting_started_with_JuMP + +A part of JuMP is `MathOptInterface` that hides the differences +behind a common interface. + +Let's try to solve some simple problem using IPOpt which is continuous +nonlinear programming solver. + +The problem is +```math +minimize \sum_{i=1}^{N}x_i^2 +``` +which is convex, so we should be able to converge to the solution +from any initial guess. + +The actual program in Julia would be +```julia +using JuMP +using Ipopt + +model = Model(Ipopt.Optimizer) +set_attrbitute(model, "output_flag", false) + +# Variables +@variable(model, x) + +@objective(model, Min, sum(x.^2)) + +optimize!(model) +``` + +One can check the final status of the optimization process via +``` +is_solved_and_feasible(model) + +termination_status(model) +``` +Here, `termination_status` returns constants from the `MathOptInterface`: +there are constants like `OPTIMAL` (global solution), +`LOCALLY_SOLVED` (local minimum), `INFEASIBLE`, etc. + +Statements for defining optimizing variables are macros: +``` +@variable(model, x) +``` +so here `x` is a symbol, and the whole thing is transformed to actual code +(probably, something like `x = variable(model, 'x')`; they do not explain it). + +To set constraints or vector variables: +``` +@variable(model, -5 <= x[i=1:42] <= 7) +``` +Also, upper and lower bounds can be set via keywords arguments to this macro. + +Interestingly, when I do +``` +typeof(x) +``` +then `x` is `VariableRef`. + + +## SciPy Optimize + +It has multiple solvers, although they have slightly different interfaces +and features: some are working with constrained optimization, some do not. + +Solution happens through a single function: +```python +from scipy import optimize + +x = minimize( + f, # Objective function + x0, # Initial guess + method="method-name", # Methods are below + args=(a, b, c, ...), # Args that are passed unfolded to f + jac=None, # Callback | True, if f returns obj and jac together + hess=None, # Callback for Hessian matrix + hessp=None, # Callback for computing Hessian-vector product + options={}, # Dictionary of options +) +``` + +The interface is general, and not all solvers (methods) use all arguments. + +Method names are case-insensitive. + +Solvers that use only the objective: +- `Nelder-Mead` +- `Powell` + +Solvers that use objective and Jacobian: +- `BFGS` Broyden-Fletcher-Goldfarb-Shanno + (can estimate Jacobian via finite differences) + +Take also Hessian or hessian product +- `Newton-CG` (Newton Conjugate Gradient) H or Hp +- `trust-NCG` Trust-region Newton-Conjugate Gradient +- `trust-Krylov` + +Only Hessian: +- `trust-exact` Trust-region Nearly Exact Algorithm: decomposes Hessian via + Cholesky factorization + +### Scipy Optimize: Constrained minimization + +Methods are: +- `trust-constr` +- `SLSQP` +- `COBYLA` +- `COBYQA` + +**Complication**: they use different interfaces to specify constraints: +`SLSQP` uses a dictionary, while others `LinearConstraint` and +`NOnlinearConstraint` instances. + +Linear constraints are written as +\[ +\begin{matrix} +c_1^\ell \\ +\dots +c_L^\ell +\end{matrix} +\leq +A x +\leq +\begin{matrix} +c_1^u \\ +\dots +c_L^u +\end{matrix} +\] diff --git a/docs/source/technotes/2025-12-10-out-and-return-args.md b/docs/source/technotes/2025-12-10-out-and-return-args.md new file mode 100644 index 000000000..586d86a16 --- /dev/null +++ b/docs/source/technotes/2025-12-10-out-and-return-args.md @@ -0,0 +1,18 @@ +# Extended Dispatch/Bridge API + +So far, we have an API like this: +``` +int call_interface_impl(implh, func_name, in_args, out_args); +``` +where `in_args` and `out_args` are arrays of input and output arguments. +Output arguments are provided by the caller and the callee writes into them. + +Now we want an extended API: +``` +int call_interface_impl(implh, func_name, in_args, out_args, return_args) +``` +where `return_args` is a semi-filled array: it contains the data types +and the number of arguments, but it is Bridge that allocates the memory +and writes into this arguments passing the ownership back to the caller +(a Converter to be precise) that converts it the return data again +from the C intermediate representation to the native language of the user. diff --git a/examples/lang_python/call_optim_convex.py b/examples/lang_python/call_optim_convex.py new file mode 100644 index 000000000..133350462 --- /dev/null +++ b/examples/lang_python/call_optim_convex.py @@ -0,0 +1,49 @@ +import argparse +import dataclasses +import sys + +import numpy as np +from openinterfaces.interfaces.optim import Optim + + +@dataclasses.dataclass +class Args: + no_artefacts: bool + + +def parse_args(argv): + p = argparse.ArgumentParser() + p.add_argument("--no-artefacts", action="store_true", help="Do not save artifacts") + args = Args(**vars(p.parse_args())) + return args + + +def main(argv=None): + if argv is None: + argv = sys.argv + args = parse_args(argv) + + x0 = np.array([1.3, 0.7, 0.8, 1.9, 1.2]) + + s = Optim("scipy_optimize") + s.set_initial_guess(x0) + s.set_objective_fn(objective_rosen_with_args) + + status, message = s.minimize() + + print(f"Status code: {status}") + print(f"Solver message: '{message}'") + print(f"Optimized value: {s.x}") + + if not args.no_artefacts: + print("Finish") + + +def objective_rosen_with_args(x, params=(0.5, 1.0)): + """The Rosenbrock function with additional arguments""" + a, b = params + return sum(a * (x[1:] - x[:-1] ** 2.0) ** 2.0 + (1 - x[:-1]) ** 2.0) + b + + +if __name__ == "__main__": + sys.exit(main()) diff --git a/include/oif/internal/bridge_api.h b/include/oif/internal/bridge_api.h index f64d848c4..611ecdf65 100644 --- a/include/oif/internal/bridge_api.h +++ b/include/oif/internal/bridge_api.h @@ -30,7 +30,8 @@ ImplInfo * load_impl(const char *impl_details, size_t version_major, size_t version_minor); int -call_impl(ImplInfo *impl_info, const char *method, OIFArgs *in_args, OIFArgs *out_args); +call_impl(ImplInfo *impl_info, const char *method, OIFArgs *in_args, OIFArgs *out_args, + OIFArgs *return_args); int unload_impl(ImplInfo *impl_info); diff --git a/include/oif/internal/dispatch.h b/include/oif/internal/dispatch.h index c7cb7e700..ab7a8ce6d 100644 --- a/include/oif/internal/dispatch.h +++ b/include/oif/internal/dispatch.h @@ -39,11 +39,13 @@ unload_interface_impl(ImplHandle implh); * @param implh Implementat handle that identifies the implementation * @param method Name of the method (function) to invoke * @param in_args Array of input arguments - * @param out_args Array of output arguments + * @param out_args Array of output (mutable) arguments + * @param return_args Array of return arguments * @return status code that signals about an error if non-zero */ int -call_interface_impl(ImplHandle implh, const char *method, OIFArgs *in_args, OIFArgs *out_args); +call_interface_impl(ImplHandle implh, const char *method, OIFArgs *in_args, OIFArgs *out_args, + OIFArgs *return_args); #ifdef __cplusplus } diff --git a/lang_c/oif_impl/bridge_c.c b/lang_c/oif_impl/bridge_c.c index 02639ed16..981f08f52 100644 --- a/lang_c/oif_impl/bridge_c.c +++ b/lang_c/oif_impl/bridge_c.c @@ -139,8 +139,10 @@ unload_impl(ImplInfo *impl_info_) } int -call_impl(ImplInfo *impl_info_, const char *method, OIFArgs *in_args, OIFArgs *out_args) +call_impl(ImplInfo *impl_info_, const char *method, OIFArgs *in_args, OIFArgs *out_args, + OIFArgs *return_args) { + (void)return_args; int result = 1; ffi_cif cif; ffi_type **arg_types = NULL; diff --git a/lang_c/oif_interfaces/src/ivp.c b/lang_c/oif_interfaces/src/ivp.c index bf5975357..942ac40a4 100644 --- a/lang_c/oif_interfaces/src/ivp.c +++ b/lang_c/oif_interfaces/src/ivp.c @@ -24,7 +24,7 @@ oif_ivp_set_rhs_fn(ImplHandle implh, oif_ivp_rhs_fn_t rhs) .arg_values = out_arg_values, }; - int status = call_interface_impl(implh, "set_rhs_fn", &in_args, &out_args); + int status = call_interface_impl(implh, "set_rhs_fn", &in_args, &out_args, NULL); return status; } @@ -48,7 +48,7 @@ oif_ivp_set_initial_value(ImplHandle implh, OIFArrayF64 *y0, double t0) .arg_values = out_arg_values, }; - int status = call_interface_impl(implh, "set_initial_value", &in_args, &out_args); + int status = call_interface_impl(implh, "set_initial_value", &in_args, &out_args, NULL); return status; } @@ -78,7 +78,7 @@ oif_ivp_set_user_data(ImplHandle implh, void *user_data) .arg_values = out_arg_values, }; - int status = call_interface_impl(implh, "set_user_data", &in_args, &out_args); + int status = call_interface_impl(implh, "set_user_data", &in_args, &out_args, NULL); return status; } @@ -101,7 +101,7 @@ oif_ivp_set_tolerances(ImplHandle implh, double rtol, double atol) .arg_values = out_arg_values, }; - int status = call_interface_impl(implh, "set_tolerances", &in_args, &out_args); + int status = call_interface_impl(implh, "set_tolerances", &in_args, &out_args, NULL); return status; } @@ -125,7 +125,7 @@ oif_ivp_integrate(ImplHandle implh, double t, OIFArrayF64 *y) .arg_values = out_arg_values, }; - int status = call_interface_impl(implh, "integrate", &in_args, &out_args); + int status = call_interface_impl(implh, "integrate", &in_args, &out_args, NULL); return status; } @@ -154,7 +154,7 @@ oif_ivp_set_integrator(ImplHandle implh, char *integrator_name, OIFConfigDict *d .arg_values = out_arg_values, }; - int status = call_interface_impl(implh, "set_integrator", &in_args, &out_args); + int status = call_interface_impl(implh, "set_integrator", &in_args, &out_args, NULL); return status; } @@ -173,6 +173,6 @@ oif_ivp_print_stats(ImplHandle implh) .arg_values = NULL, }; - int status = call_interface_impl(implh, "print_stats", &in_args, &out_args); + int status = call_interface_impl(implh, "print_stats", &in_args, &out_args, NULL); return status; } diff --git a/lang_c/oif_interfaces/src/linsolve.c b/lang_c/oif_interfaces/src/linsolve.c index eea129717..b97337027 100644 --- a/lang_c/oif_interfaces/src/linsolve.c +++ b/lang_c/oif_interfaces/src/linsolve.c @@ -21,7 +21,7 @@ oif_solve_linear_system(ImplHandle implh, OIFArrayF64 *A, OIFArrayF64 *b, OIFArr .arg_values = out_arg_values, }; - int status = call_interface_impl(implh, "solve_lin", &in_args, &out_args); + int status = call_interface_impl(implh, "solve_lin", &in_args, &out_args, NULL); return status; } diff --git a/lang_c/oif_interfaces/src/qeq.c b/lang_c/oif_interfaces/src/qeq.c index f6f6e5578..410331eeb 100644 --- a/lang_c/oif_interfaces/src/qeq.c +++ b/lang_c/oif_interfaces/src/qeq.c @@ -21,7 +21,7 @@ oif_solve_qeq(ImplHandle implh, double a, double b, double c, OIFArrayF64 *roots .arg_values = out_arg_values, }; - int status = call_interface_impl(implh, "solve_qeq", &in_args, &out_args); + int status = call_interface_impl(implh, "solve_qeq", &in_args, &out_args, NULL); return status; } diff --git a/lang_julia/OpenInterfaces/src/OpenInterfaces.jl b/lang_julia/OpenInterfaces/src/OpenInterfaces.jl index f7faca3cf..073aa1b00 100644 --- a/lang_julia/OpenInterfaces/src/OpenInterfaces.jl +++ b/lang_julia/OpenInterfaces/src/OpenInterfaces.jl @@ -160,6 +160,7 @@ function call_impl( func_name::String, in_user_args::Tuple{Vararg{Any}}, out_user_args::Tuple{Vararg{Any}}, + return_user_args::Tuple{Vararg{Any}} = (), )::Int in_num_args = length(in_user_args) out_num_args = length(out_user_args) @@ -299,6 +300,10 @@ function call_impl( in_args = Ref(OIFArgs(in_num_args, pointer(in_arg_types), pointer(in_arg_values))) out_args = Ref(OIFArgs(out_num_args, pointer(out_arg_types), pointer(out_arg_values))) + if (isempty(return_user_args)) + return_args = C_NULL + end + result = GC.@preserve in_arg_types in_arg_values out_arg_types out_arg_values temp_refs begin @ccall $(call_interface_impl_fn[])( @@ -306,6 +311,7 @@ function call_impl( func_name::Cstring, in_args::Ptr{OIFArgs}, out_args::Ptr{OIFArgs}, + return_args::Ptr{OIFArgs}, )::Int end diff --git a/lang_julia/bridge_julia.c b/lang_julia/bridge_julia.c index 92e1188e1..ac7594485 100644 --- a/lang_julia/bridge_julia.c +++ b/lang_julia/bridge_julia.c @@ -480,8 +480,10 @@ unload_impl(ImplInfo *impl_info_) } int -call_impl(ImplInfo *impl_info_, const char *method, OIFArgs *in_args, OIFArgs *out_args) +call_impl(ImplInfo *impl_info_, const char *method, OIFArgs *in_args, OIFArgs *out_args, + OIFArgs *return_args) { + (void)return_args; int result = -1; if (impl_info_->dh != OIF_LANG_JULIA) { diff --git a/lang_python/oif/openinterfaces/core.py b/lang_python/oif/openinterfaces/core.py index 44b0fd8b2..9433cb19d 100644 --- a/lang_python/oif/openinterfaces/core.py +++ b/lang_python/oif/openinterfaces/core.py @@ -105,6 +105,8 @@ else: _lib_dispatch = ctypes.PyDLL(_lib_dispatch_name) +_liboif_common_util = ctypes.PyDLL("liboif_common_util.so") + elapsed = 0.0 elapsed_call = 0.0 @@ -249,6 +251,11 @@ def wrapper(*arg_values): def make_oif_user_data(data: object) -> OIFUserData: + """Make OIFUserData structure for Python object `data`. + + *IMPORTANT* Keep the reference to `data`, otherwise it can be + garbage collected, and then using it becomes a segfault situation. + """ return OIFUserData(OIF_LANG_PYTHON, None, None, ctypes.c_void_p(id(data))) @@ -297,7 +304,7 @@ def __init__(self, implh, interface, impl): ], ) - def call(self, method, user_args, out_user_args): + def call(self, method, user_args, out_user_args, return_arg_types=()): num_args = len(user_args) arg_types = [] arg_values = [] @@ -419,17 +426,59 @@ def call(self, method, user_args, out_user_args): ) out_packed = OIFArgs(num_out_args, out_arg_types_ctypes, out_arg_values_ctypes) + # Process return arguments, that is the arguments + # that the callee should allocate and return to us. + return_args_packed = None + if return_arg_types: + num_return_args = len(return_arg_types) + return_arg_types_ctypes = ctypes.cast( + (ctypes.c_int * len(return_arg_types))(*return_arg_types), + ctypes.POINTER(OIFArgType), + ) + return_arg_values = [None] * len(return_arg_types) + return_arg_values_ctypes = ctypes.cast( + (ctypes.c_void_p * len(return_arg_types))(*return_arg_values), + ctypes.POINTER(ctypes.c_void_p), + ) + return_args_packed = ctypes.pointer( + OIFArgs( + num_return_args, return_arg_types_ctypes, return_arg_values_ctypes + ) + ) + status = self._call_interface_impl( self.implh, method.encode(), ctypes.byref(in_args_packed), ctypes.byref(out_packed), + return_args_packed, ) if status != 0: raise RuntimeError(f"Error occurred while executing method '{method}'") - return 0 + return_args = [] + if return_args_packed is not None: + oif_util_free_ = _wrap_c_function( + _liboif_common_util, "oif_util_free_", None, [ctypes.c_void_p] + ) + + for i, arg_type in enumerate(return_arg_types): + void_pointer = ctypes.c_void_p( + return_args_packed.contents.arg_values[i] + ) + if arg_type == OIF_TYPE_INT: + py_var = ctypes.cast( + void_pointer, ctypes.POINTER(ctypes.c_int) + ).contents.value + elif arg_type == OIF_TYPE_STRING: + py_var = ctypes.string_at(void_pointer).decode("utf-8") + else: + raise RuntimeError("Cannot convert return argument") + return_args.append(py_var) + oif_util_free_(void_pointer.value) + + return return_args if return_args else status def load_impl(interface: str, impl: str, major: UInt, minor: UInt): diff --git a/lang_python/oif_impl/bridge_python.c b/lang_python/oif_impl/bridge_python.c index 7f0ecf815..46d41bd08 100644 --- a/lang_python/oif_impl/bridge_python.c +++ b/lang_python/oif_impl/bridge_python.c @@ -30,6 +30,8 @@ static int IMPL_COUNTER = 0; static bool is_python_initialized_by_us = false; +static bool NUMPY_IS_INITIALIZED_ = false; + static char prefix_[] = "bridge_python"; PyObject * @@ -196,35 +198,23 @@ init_python_() return 0; } -ImplInfo * -load_impl(const char *impl_details, size_t version_major, size_t version_minor) +static void * +init_numpy_() { - PyObject *pFileName, *pModule; - PyObject *pClass, *pInstance; - PyObject *pFunc; - PyObject *pInitArgs; - PyObject *pArgs; - PyObject *pValue; + PyObject *pFileName = NULL; + PyObject *pModule = NULL; + PyObject *pClass = NULL; + PyObject *pInstance = NULL; + PyObject *pFunc = NULL; + PyObject *pInitArgs = NULL; + PyObject *pArgs = NULL; + PyObject *pValue = NULL; - ImplInfo *result = NULL; int status; int nbytes_required; // Number of required bytes in snprintf. int nbytes_written; // Number of written bytes in snprintf. char *cmd; // Command formatted in snprintf. - - (void)version_major; - (void)version_minor; - if (Py_IsInitialized()) { - fprintf(stderr, "[%s] Backend is already initialized\n", prefix_); - } - else { - int status; - status = init_python_(); - if (status != 0) { - return NULL; - } - } - + // // We need to `dlopen` the Python library, otherwise, // NumPy initialization fails. // Details: @@ -299,12 +289,18 @@ load_impl(const char *impl_details, size_t version_major, size_t version_minor) } status = PyRun_SimpleString(cmd); free(cmd); + cmd = NULL; if (status < 0) { logerr(prefix_, "An error occurred when initializating Python"); goto cleanup; } - import_array2("Failed to initialize NumPy C API", NULL); + status = _import_array(); + if (status != 0) { + PyErr_Print(); + logerr(prefix_, "Failed to initialize NumPy C API"); + goto cleanup; + } const char *cmd_numpy_fmt = "import numpy; " @@ -319,11 +315,64 @@ load_impl(const char *impl_details, size_t version_major, size_t version_minor) } status = PyRun_SimpleString(cmd); free(cmd); + cmd = NULL; if (status < 0) { logerr(prefix_, "An error occurred when initializating Python"); goto cleanup; } + NUMPY_IS_INITIALIZED_ = true; + +cleanup: + Py_XDECREF(pValue); + Py_XDECREF(pArgs); + Py_XDECREF(pInitArgs); + Py_XDECREF(pFunc); + Py_XDECREF(pInstance); + Py_XDECREF(pClass); + // Somehow, the reference to sysconfig must be preserved. + // Otherwise, there will be a crash when loading implementations. + /* Py_XDECREF(pModule); */ + Py_XDECREF(pFileName); + + if (cmd != NULL) { + free(cmd); + } + return NULL; +} + +ImplInfo * +load_impl(const char *impl_details, size_t version_major, size_t version_minor) +{ + PyObject *pFileName = NULL; + PyObject *pModule = NULL; + PyObject *pClass = NULL; + PyObject *pInstance = NULL; + PyObject *pFunc = NULL; + PyObject *pInitArgs = NULL; + PyObject *pArgs = NULL; + PyObject *pValue = NULL; + + ImplInfo *result = NULL; + int status; + + (void)version_major; + (void)version_minor; + if (Py_IsInitialized()) { + fprintf(stderr, "[%s] Backend is already initialized\n", prefix_); + } + else { + int status; + status = init_python_(); + if (status != 0) { + return NULL; + } + } + + if (!NUMPY_IS_INITIALIZED_) { + init_numpy_(); + } + char moduleName[512] = "\0"; char className[512] = "\0"; size_t i; @@ -378,7 +427,6 @@ load_impl(const char *impl_details, size_t version_major, size_t version_minor) logerr(prefix_, "Failed to instantiate class %s\n", className); goto cleanup; } - /* Py_INCREF(pInstance); */ PythonImplInfo *impl_info = oif_util_malloc(sizeof(*impl_info)); if (impl_info == NULL) { @@ -409,7 +457,8 @@ load_impl(const char *impl_details, size_t version_major, size_t version_minor) } int -call_impl(ImplInfo *impl_info, const char *method, OIFArgs *in_args, OIFArgs *out_args) +call_impl(ImplInfo *impl_info, const char *method, OIFArgs *in_args, OIFArgs *out_args, + OIFArgs *return_args) { int result = 1; if (impl_info->dh != OIF_LANG_PYTHON) { @@ -555,6 +604,10 @@ call_impl(ImplInfo *impl_info, const char *method, OIFArgs *in_args, OIFArgs *ou else if (out_args->arg_types[i] == OIF_TYPE_ARRAY_F64) { pValue = get_numpy_array_from_oif_array_f64(out_args->arg_values[i]); } + else if (out_args->arg_types[i] == OIF_TYPE_STRING) { + pValue = PyArray_SimpleNewFromData(1, (intptr_t[1]){1000}, NPY_UINT8, + *(char **)out_args->arg_values[i]); + } else { pValue = NULL; } @@ -569,14 +622,75 @@ call_impl(ImplInfo *impl_info, const char *method, OIFArgs *in_args, OIFArgs *ou // Invoke function. pValue = PyObject_CallObject(pFunc, pArgs); Py_DECREF(pArgs); - if (pValue != NULL) { - Py_DECREF(pValue); - } - else { - Py_DECREF(pFunc); + if (pValue == NULL) { PyErr_Print(); - fprintf(stderr, "[%s] Call failed\n", prefix_); - return 2; + logerr(prefix_, "Call to the method '%s' has failed\n", method); + goto cleanup; + } + + if (return_args != NULL) { + if (return_args->num_args > 1 && !PyTuple_Check(pValue)) { + logerr(prefix_, + "Expected %d return arguments, however the method '%s' " + "has not returned a tuple", + return_args->num_args, method); + goto cleanup; + } + if (return_args->num_args != 1 && return_args->num_args != PyTuple_Size(pValue)) { + logerr(prefix_, + "Expected %d returned arguments, " + "got %d instead in the call to the method '%s'", + return_args->num_args, PyTuple_Size(pValue), method); + goto cleanup; + } + + for (size_t i = 0; i < return_args->num_args; ++i) { + PyObject *val = PyTuple_GetItem(pValue, i); + switch (return_args->arg_types[i]) { + case OIF_TYPE_I32: + if (!PyLong_Check(val)) { + logerr(prefix_, + "Expected Python integer object, but did not get it"); + goto cleanup; + } + + long tmp = PyLong_AsLong(val); + if (tmp == -1 && PyErr_Occurred()) { + logerr(prefix_, "Could not convert Python object to C long value"); + goto cleanup; + } + if (tmp >= INT32_MIN && tmp <= INT32_MAX) { + return_args->arg_values[i] = oif_util_malloc(sizeof(int32_t)); + *(int32_t *)return_args->arg_values[i] = tmp; + } + else { + logerr(prefix_, "Return value is outside of the range of int32"); + goto cleanup; + } + break; + case OIF_TYPE_STRING: + if (!PyUnicode_Check(val)) { + logerr(prefix_, + "Expected Unicode string object, but did not get it"); + goto cleanup; + } + // Quote from docs for `PyUnicode_AsUTF8`: + // The caller is not responsible for deallocating the buffer. + const char *s = PyUnicode_AsUTF8(val); + if (s == NULL) { + logerr(prefix_, + "Expected Unicode string object, but did not get it"); + } + size_t length = PyUnicode_GET_LENGTH(val) + 1; + return_args->arg_values[i] = oif_util_malloc(sizeof(char) * length); + snprintf(return_args->arg_values[i], length, "%s", s); + break; + default: + logerr(prefix_, "Return type with id '%d' is not supported yet", + return_args->arg_types[i]); + goto cleanup; + } + } } } else { diff --git a/lang_python/oif_impl/openinterfaces/_impl/interfaces/optim.py b/lang_python/oif_impl/openinterfaces/_impl/interfaces/optim.py new file mode 100644 index 000000000..d220077be --- /dev/null +++ b/lang_python/oif_impl/openinterfaces/_impl/interfaces/optim.py @@ -0,0 +1,36 @@ +import abc +from typing import Callable, Union + +import numpy as np + + +class OptimInterface(abc.ABC): + @abc.abstractmethod + def set_initial_guess(self, x0: np.ndarray) -> Union[int, None]: + """Set initial guess for minimization.""" + + @abc.abstractmethod + def set_objective_fn(self, objective_fn: Callable) -> Union[int, None]: + """Specify objective function f.""" + + # @abc.abstractmethod + # def set_tolerances(self, rtol: float, atol: float) -> Union[int, None]: + # """Specify relative and absolute tolerances, respectively.""" + + # @abc.abstractmethod + # def set_user_data(self, user_data: object) -> Union[int, None]: + # """Specify additional data that will be used for right-hand side function.""" + + # @abc.abstractmethod + # def set_integrator( + # self, integrator_name: str, integrator_params: Dict + # ) -> Union[int, None]: + # """Set integrator, if the name is recognizable.""" + + # @abc.abstractmethod + # def integrate(self, t: float, y: np.ndarray) -> Union[int, None]: + # """Integrate to time `t` and write solution to `y`.""" + + # @abc.abstractmethod + # def print_stats(self): + # """Print integration statistics.""" diff --git a/lang_python/oif_impl/openinterfaces/_impl/optim/scipy_optimize/scipy_optimize.conf b/lang_python/oif_impl/openinterfaces/_impl/optim/scipy_optimize/scipy_optimize.conf new file mode 100644 index 000000000..dd2eebf0a --- /dev/null +++ b/lang_python/oif_impl/openinterfaces/_impl/optim/scipy_optimize/scipy_optimize.conf @@ -0,0 +1,2 @@ +python +openinterfaces._impl.optim.scipy_optimize.scipy_optimize ScipyOptimize 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 new file mode 100644 index 000000000..a8d6d422f --- /dev/null +++ b/lang_python/oif_impl/openinterfaces/_impl/optim/scipy_optimize/scipy_optimize.py @@ -0,0 +1,106 @@ +from typing import Callable + +import numpy as np +from openinterfaces._impl.interfaces.optim import OptimInterface +from scipy import optimize + +_prefix = "scipy_optimize" + + +class ScipyOptimize(OptimInterface): + """ + Solve optimization problem minimize_x f(x) + + """ + + def __init__(self) -> None: + self.N = 0 # Problem dimension. + self.x0: np.ndarray | None = None + self.objective_fn: Callable | None = None + self.user_data: object | None = None + self.with_user_data = False + self.method_name: str | None = None + self.method_params: str | None = None + + def set_initial_guess(self, x0: np.ndarray): + self.x0 = x0 + self.N = len(x0) + + def set_user_data(self, user_data): + self.user_data = user_data + self.with_user_data = True + + def set_objective_fn(self, objective_fn): + self.objective_fn = objective_fn + + x = np.random.random(size=(self.N,)) + msg = "Wrong signature for the objective function: " + msg += "expected return value must be of `float64` type" + + if self.with_user_data: + assert type(self.objective_fn(x, self.user_data)) in [ + float, + np.float64, + ], msg + else: + assert type(self.objective_fn(x)) in [float, np.float64], msg + + def set_method(self, method_name, method_params): + available_options = optimize.show_options( + "minimize", method=method_name, disp=False + ) + + for k in method_params.keys(): + candidates = [f"\n{k} :", f"\n{k},"] + found = False + for c in candidates: + if c in available_options: + found = True + break + + if not found: + raise ValueError( + f"Method option '{k}' is not known for the method '{method_name}'" + ) + + self.method_name = method_name + self.method_params = method_params + + def minimize(self, out_x): + if self.x0 is None: + raise RuntimeError("Method `set_initial_guess` must be called first") + + if self.objective_fn is None: + raise RuntimeError("Method `set_objective_fn` must be called first") + + if len(self.x0) != len(out_x): + raise ValueError( + "Shapes of the output array and the initial-guess array differ" + ) + + if self.with_user_data: + result = optimize.minimize( + self.objective_fn, + self.x0, + args=self.user_data, + method=self.method_name, + options=self.method_params, + ) + else: + result = optimize.minimize( + self.objective_fn, + self.x0, + method=self.method_name, + options=self.method_params, + ) + + out_x[:] = result.x + return (result.status, result.message) + + # def print_stats(self): + # print("WARNING: `scipy_ode` does not provide statistics") + + # def _rhs_fn_wrapper(self, t, y): + # """Callback that satisfies signature expected by Open Interfaces.""" + # self.rhs(t, y, self.ydot, self.user_data) + # return self.ydot diff --git a/lang_python/oif_interfaces/openinterfaces/interfaces/optim.py b/lang_python/oif_interfaces/openinterfaces/interfaces/optim.py new file mode 100644 index 000000000..f8b9b59bd --- /dev/null +++ b/lang_python/oif_interfaces/openinterfaces/interfaces/optim.py @@ -0,0 +1,144 @@ +r"""This module defines the interface for solving minimization problems: + +.. math:: + minimize_x f(x) + +where :math:`f : \mathbb R^n \to \mathbb R`. +""" + +from collections.abc import Callable +from typing import TypeAlias + +import numpy as np +from openinterfaces.core import ( + OIF_TYPE_ARRAY_F64, + OIF_TYPE_F64, + OIF_TYPE_INT, + OIF_TYPE_STRING, + OIF_USER_DATA, + OIFPyBinding, + OIFUserData, + load_impl, + make_oif_callback, + make_oif_user_data, + unload_impl, +) + +ObjectiveFn: TypeAlias = Callable[[np.ndarray, object], int] +"""Signature of the objective function :math:`f(, y)`. + +!!!!!!!!! The function accepts four arguments: +!!!!!!!!! - `t`: current time, +!!!!!!!!! - `y`: state vector at time :math:`t`, +!!!!!!!!! - `ydot`: output array to which the result of ??? +!!!!!!!!! - `user_data`: additional context (user-defined data) that +!!!!!!!!! must be passed to the function (e.g., parameters of the system). + +""" + + +class Optim: + r"""Interface for solving optimization (minimization) problems. + + This class serves as a gateway to the implementations of the + solvers for optimization problems. + + Parameters + ---------- + impl : str + Name of the desired implementation. + + Examples + -------- + + Let's solve the following convex optimization problem: + + .. math:: + minimize \sum_{i = 1}^N x_i^2 + + where the solution is :math:`[0, ..., 0]` as the problem is convex. + + First, import the necessary modules: + >>> import numpy as np + >>> from oif.interfaces.optim import Optim + + Define the objective function: + + >>> def objective_fn(x): + ... return np.sum(x**2) + + + Create an instance of the optim solver using the implementation "scipy_optimize", + which is an adapter to the `scipy.optimize` Python package: + + >>> s = Optim("scipy_optimize") + + We set the initial value, the right-hand side function, and the tolerance: + + >>> s.set_initial_guess([2.718, 3.142]) + >>> s.set_objective_fn(objective_fn) + + Now we solve the minimization problem and print the return status and message: + + >>> status, message = s.minimize() + + We can print the resultant minimizer by retrieving it from the solver: + + >>> print(f"Minimizer is {s.x}") + + """ + + def __init__(self, impl: str): + self._binding: OIFPyBinding = load_impl("optim", impl, 1, 0) + self.x0: np.ndarray + """Current value of the state vector.""" + self._N: int = 0 + self.x: np.ndarray + self.user_data: object + self.oif_user_data: OIFUserData + + def set_initial_guess(self, x0: np.ndarray): + """Set initial guess for the optimization problem""" + self.x0 = x0 + self._N = len(self.x0) + self.x = np.empty((self._N,)) + + # It is not clear to me whether we should pass `self.x0` + # and therefore keep the reference to `x0` to avoid crashes + # or an implementation is supposed to copy it. + self._binding.call("set_initial_guess", (self.x0,), ()) + + def set_objective_fn(self, objective_fn: ObjectiveFn): + self.wrapper = make_oif_callback( + objective_fn, (OIF_TYPE_ARRAY_F64, OIF_USER_DATA), OIF_TYPE_F64 + ) + self._binding.call("set_objective_fn", (self.wrapper,), ()) + + def set_user_data(self, user_data: object): + """Specify additional data that will be used for right-hand side function.""" + self.user_data = user_data + self.oif_user_data = make_oif_user_data(user_data) + self._binding.call("set_user_data", (self.oif_user_data,), ()) + + def set_method(self, method_name: str, method_params: dict = {}): + """Set integrator, if the name is recognizable.""" + self._binding.call("set_method", (method_name, method_params), ()) + + def minimize(self): + """Integrate to time `t` and write solution to `y`.""" + status, message = self._binding.call( + "minimize", + (), + (self.x,), + (OIF_TYPE_INT, OIF_TYPE_STRING), + ) + + return status, message + + # def print_stats(self): + # """Print integration statistics.""" + # self._binding.call("print_stats", (), ()) + + def __del__(self): + if hasattr(self, "_binding"): + unload_impl(self._binding) diff --git a/tests/lang_python/test_optim.py b/tests/lang_python/test_optim.py new file mode 100644 index 000000000..5abe70ab5 --- /dev/null +++ b/tests/lang_python/test_optim.py @@ -0,0 +1,95 @@ +import numpy as np +import pytest +from openinterfaces.interfaces.optim import Optim + + +# ----------------------------------------------------------------------------- +# Problems +def convex_objective_fn(x): + return np.sum(x**2) + + +def convex_objective_with_args_fn(x, args): + return np.sum((x - args) ** 2) + + +def rosenbrock_objective_fn(x): + """The Rosenbrock function with additional arguments""" + a, b = (0.5, 1.0) + return sum(a * (x[1:] - x[:-1] ** 2.0) ** 2.0 + (1 - x[:-1]) ** 2.0) + b + + +@pytest.fixture(params=["scipy_optimize"]) +def s(request): + return Optim(request.param) + + +# ----------------------------------------------------------------------------- +# Tests +def test__simple_convex_problem__converges(s): + x0 = np.array([0.5, 0.6, 0.7]) + + s.set_initial_guess(x0) + s.set_objective_fn(convex_objective_fn) + + status, message = s.minimize() + x = s.x + + assert status == 0 + assert isinstance(message, str) and len(message) > 0 + assert len(x) == len(x0) + assert np.all(np.abs(x) < 1e-6) + + +def test__rosenbrock_problem__converges(s): + x0 = np.array([3.14, 2.72, 42.0, 9.81, 8.31]) + + s.set_initial_guess(x0) + s.set_objective_fn(rosenbrock_objective_fn) + + status, message = s.minimize() + x = s.x + + assert status == 0 + assert len(x) == len(x0) + assert np.all(np.abs(x - 1.0) < 1e-5) # The solution is [1, 1, ..., 1]. + + +def test__parameterized_convex_problem__converges(s): + x0 = np.array([0.5, 0.6, 0.7]) + user_data = np.array([2.0, 3.0, -1.0]) + + s.set_initial_guess(x0) + s.set_user_data(user_data) + s.set_objective_fn(convex_objective_with_args_fn) + + status, message = s.minimize() + x = s.x + + assert status == 0 + assert len(x) == len(x0) + 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]) + + 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) + + +def test__set_method_with__wrong_method_params__are_not_accepted(s): + with pytest.raises(RuntimeError): + s.set_method("nelder-mead", {"wrong_param": 42})