diff --git a/.github/workflows/build-python-wheels.yml b/.github/workflows/build-python-wheels.yml index 1dd7eb5cc..99f65bce9 100644 --- a/.github/workflows/build-python-wheels.yml +++ b/.github/workflows/build-python-wheels.yml @@ -148,7 +148,6 @@ jobs: CMAKE_Fortran_COMPILER=C:/mingw64/bin/gfortran.exe CIBW_CONFIG_SETTINGS_WINDOWS: > cmake.define.METIS_LIBRARY=dependencies/bin/libmetis.dll - cmake.define.BLAS_LIBRARIES="dependencies/lib/libblas.a;dependencies/lib/libcblas.a" cmake.define.AUXILIARY_LIBRARIES=dependencies/lib/libstdc++.a - name: Store artifacts @@ -287,4 +286,4 @@ jobs: - name: Run example working-directory: ${{github.workspace}}/interfaces/Python/example - run: python example_hs015.py \ No newline at end of file + run: python example_hs015.py diff --git a/.github/workflows/matlab-test-example.yml b/.github/workflows/matlab-test-example.yml new file mode 100644 index 000000000..66cc3f87b --- /dev/null +++ b/.github/workflows/matlab-test-example.yml @@ -0,0 +1,104 @@ +name: Test the Matlab interface on the example + +on: + push: + branches: [ "main" ] + paths-ignore: + - '*.md' + - 'LICENSE' + - '*.cff' + - '*.yml' + - '*.yaml' + - 'docs/**' + pull_request: + branches: [ "main" ] + paths-ignore: + - '*.md' + - 'LICENSE' + - '*.cff' + - '*.yml' + - '*.yaml' + - 'docs/**' + +env: + BUILD_TYPE: Debug + +jobs: + build: + runs-on: ${{ matrix.os }} + strategy: + fail-fast: false + matrix: + include: + - os: windows-latest + addpath: install/bin + - os: macos-latest + addpath: install/lib + + steps: + - uses: actions/checkout@v4 + + - name: Set up MATLAB + uses: matlab-actions/setup-matlab@v2 + + - name: Install Fortran compiler (macOS) + if: startsWith(matrix.os, 'macos') + uses: fortran-lang/setup-fortran@main + with: + compiler: 'gcc' + version: '14' + + - name: Setup MSYS2 (Windows) + if: startsWith(matrix.os, 'windows') + uses: msys2/setup-msys2@v2 + with: + path-type: inherit + msystem: MINGW64 + update: true + install: | + mingw-w64-x86_64-gcc-fortran + mingw-w64-x86_64-cmake + mingw-w64-x86_64-make + + - name: Download dependencies + run: bash dependencies/scripts/download_dependencies.sh + + - name: Configure (macOS) + if: startsWith(matrix.os, 'macos') + run: | + cmake -B ${{github.workspace}}/build -DCMAKE_BUILD_TYPE=${{env.BUILD_TYPE}} \ + -DBQPD=${{github.workspace}}/dependencies/lib/libbqpd.a \ + -DLAPACK_LIBRARIES=${{github.workspace}}/dependencies/lib/liblapack.a \ + -DBLAS_LIBRARIES="${{github.workspace}}/dependencies/lib/libblas.a;${{github.workspace}}/dependencies/lib/libcblas.a" + + - name: Configure (Windows) + if: startsWith(matrix.os, 'windows') + shell: cmd + run: | + cmake -G "MinGW Makefiles" -B "${{github.workspace}}/build" ^ + -DBQPD=${{github.workspace}}/dependencies/lib/libbqpd.a ^ + -DLAPACK_LIBRARIES=${{github.workspace}}/dependencies/lib/liblapack.a ^ + -DBLAS_LIBRARIES="${{github.workspace}}/dependencies/lib/libblas.a;${{github.workspace}}/dependencies/lib/libcblas.a" + + - name: Build + run: cmake --build ${{github.workspace}}/build --config ${{env.BUILD_TYPE}} --parallel + + - name: Build unomex + run: cmake --build ${{github.workspace}}/build --target unomex --config ${{env.BUILD_TYPE}} --parallel + + - name: Install unomex + run: cmake --install ${{github.workspace}}/build --prefix ${{github.workspace}}/install + + - name: Run example (high-level interface) + uses: matlab-actions/run-command@v2 + with: + command: | + addpath('${{ matrix.addpath }}'); + run('interfaces/Matlab/example/example_uno.m'); + + - name: Run example (low-level interface) + uses: matlab-actions/run-command@v2 + with: + command: | + addpath('${{ matrix.addpath }}'); + run('interfaces/Matlab/example/example_hs015.m'); \ No newline at end of file diff --git a/CMakeLists.txt b/CMakeLists.txt index c98e66635..81b42d6b5 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -370,18 +370,70 @@ if(SKBUILD) target_compile_definitions(unopy PRIVATE PYBIND11_DETAILED_ERROR_MESSAGES) endif() - if (CMAKE_CXX_COMPILER_ID STREQUAL "GNU") - # link-time optimization - target_compile_options(unopy PRIVATE $<$:-flto=auto>) - target_link_options(unopy PRIVATE -flto=auto) - endif() - install(TARGETS unopy LIBRARY DESTINATION unopy # Linux/macOS: unopy/unopy.so RUNTIME DESTINATION unopy # Windows: unopy/unopy.pyd ) endif() +############################ +# optional Matlab bindings # +############################ + +find_package(Matlab 23.2 COMPONENTS MAIN_PROGRAM MEX_COMPILER) # >=23.2 = R2023b +if(Matlab_FOUND) + message(STATUS "Found Matlab ${Matlab_VERSION}") + + file(GLOB MATLAB_SOURCE_FILES + interfaces/Matlab/cpp_classes/*.cpp + interfaces/Matlab/unomex/*.cpp + ) + + # libut matlab library + # cf. https://undocumentedmatlab.com/articles/mex-ctrl-c-interrupt#:~:text=The%20relevant%20functions%20seem%20to,utLongjmpIfInterruptPending(?) + get_filename_component(Matlab_MX_LIBRARY_PATH ${Matlab_MX_LIBRARY} DIRECTORY ) + find_library(MXLIBUT NAMES libut ut + PATHS ${Matlab_MX_LIBRARY_PATH}) + if(MXLIBUT) + add_definitions("-D HAS_MXLIBUT") + endif() + + # compile the Unomex source files once + add_library(compiled_unomex_files EXCLUDE_FROM_ALL OBJECT ${MATLAB_SOURCE_FILES}) + target_link_libraries(compiled_unomex_files PUBLIC uno_dependencies) + set_property(TARGET compiled_unomex_files PROPERTY POSITION_INDEPENDENT_CODE ON) + target_include_directories(compiled_unomex_files PUBLIC ${Matlab_INCLUDE_DIRS}) + + # uno_optimize + matlab_add_mex(NAME unomex-optimize SRC $ interfaces/Matlab/uno_optimize.cpp + EXCLUDE_FROM_ALL + OUTPUT_NAME uno_optimize + LINK_TO ${DEFAULT_UNO_LIB} uno_dependencies ${MXLIBUT}) + + # uno_options + matlab_add_mex(NAME unomex-options SRC $ interfaces/Matlab/uno_options.cpp + EXCLUDE_FROM_ALL + OUTPUT_NAME uno_options + LINK_TO ${DEFAULT_UNO_LIB} uno_dependencies ${MXLIBUT}) + target_include_directories(unomex-options PUBLIC ${directory}) + + # logical target to compile both mex + add_custom_target(unomex DEPENDS unomex-optimize unomex-options) + + # cf. https://github.com/cvanaret/Uno/pull/360#issuecomment-3431367949 + if(APPLE) + set(MEX_RPATH "@executable_path/../extern/bin/maci64;@executable_path/../extern/bin/maca64;@executable_path/../../extern/bin/maci64;@executable_path/../../extern/bin/maca64") + set_target_properties(unomex-optimize PROPERTIES BUILD_RPATH "${MEX_RPATH};${CMAKE_BUILD_RPATH}" INSTALL_RPATH "${MEX_RPATH};${CMAKE_INSTALL_RPATH}") + set_target_properties(unomex-options PROPERTIES BUILD_RPATH "${MEX_RPATH};${CMAKE_BUILD_RPATH}" INSTALL_RPATH "${MEX_RPATH};${CMAKE_INSTALL_RPATH}") + endif() + + # install + set(MATLAB_INSTALL_FILES "interfaces/Matlab/uno.m;interfaces/Matlab/uno_optimize.m;interfaces/Matlab/uno_options.m") + install(TARGETS unomex-optimize LIBRARY DESTINATION ${CMAKE_INSTALL_BINDIR}) + install(TARGETS unomex-options LIBRARY DESTINATION ${CMAKE_INSTALL_BINDIR}) + install(FILES ${MATLAB_INSTALL_FILES} DESTINATION $,bin,lib>) +endif() + # ====================== # GoogleTest unit tests # ====================== @@ -434,4 +486,4 @@ install(TARGETS ${INSTALL_TARGETS} ARCHIVE DESTINATION ${CMAKE_INSTALL_LIBDIR} # static libraries LIBRARY DESTINATION ${CMAKE_INSTALL_LIBDIR} # shared libraries PUBLIC_HEADER DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/uno # headers -) \ No newline at end of file +) diff --git a/interfaces/C/Uno_C_API.cpp b/interfaces/C/Uno_C_API.cpp index b710f2a09..f9f47eeec 100644 --- a/interfaces/C/Uno_C_API.cpp +++ b/interfaces/C/Uno_C_API.cpp @@ -3,8 +3,6 @@ #include #include -#include -#include #include "Uno_C_API.h" #include "../UserModel.hpp" #include "Uno.hpp" @@ -379,75 +377,29 @@ class CUserCallbacks: public UserCallbacks { void* user_data; }; -// std::streambuf wrapper around uno_logger_stream_callback -class CStreamBuffer : public std::streambuf { -public: - explicit CStreamBuffer(uno_logger_stream_callback logger_stream_callback, void* user_data, std::size_t buffer_size) : - logger_stream_callback(logger_stream_callback), user_data(user_data) { - // allocate output buffer and set stream buffer pointer - this->buffer = new char[buffer_size]; - this->setp(this->buffer, this->buffer + buffer_size - 1); - } - ~CStreamBuffer() override { - // flush remaining data and release buffer memory - this->sync(); - delete[] this->buffer; - } +class CStreamCallback : public UserStreamCallback { +public: + CStreamCallback(uno_logger_stream_callback logger_stream_callback, void* user_data) : + UserStreamCallback(), logger_stream_callback(logger_stream_callback), user_data(user_data) { } + ~CStreamCallback() override { } -protected: - // called on buffer overflow - int overflow(int character) override { - if (character != EOF) { - // insert the character into the buffer - *this->pptr() = traits_type::to_char_type(character); - this->pbump(1); + int32_t operator()(const char* buf, int32_t len) const override { + if (this->logger_stream_callback) { + return this->logger_stream_callback(buf, len, this->user_data); + } + else { + return 0; } - // return EOF for error - return (this->flush_buffer() == 0) ? character : EOF; - } - - int sync() override { - return this->flush_buffer(); } private: uno_logger_stream_callback logger_stream_callback; void* user_data; char* buffer; - - // flush buffer to the logger callback - int flush_buffer() { - // check for invalid stream callback - if (!this->logger_stream_callback) { - return -1; - } - std::ptrdiff_t current_used_buffer_size = this->pptr() - this->pbase(); - if (current_used_buffer_size > 0) { - // call user logger callback - const uno_int callback_result = this->logger_stream_callback(this->pbase(), static_cast(current_used_buffer_size), - this->user_data); - if (callback_result != static_cast(current_used_buffer_size)) { - return -1; - } - // move buffer pointer - this->pbump(static_cast(-current_used_buffer_size)); - } - return 0; - } -}; - -// std::ostream wrapper around uno_logger_stream_callback -class COStream : public std::ostream { -public: - COStream(uno_logger_stream_callback logger_stream_callback, void* user_data, std::size_t buffer_size = 1024) : // 1024 default buffer size - std::ostream(&this->buffer), buffer(logger_stream_callback, user_data, buffer_size) { } - -private: - // internal stream buffer that sends output to the uno_logger_stream_callback - CStreamBuffer buffer; }; -COStream* c_ostream = nullptr; +CStreamCallback* c_stream_callback = nullptr; +UserOStream* ostream = nullptr; struct Solver { Uno* solver; @@ -939,15 +891,19 @@ bool uno_set_solver_callbacks(void* solver, uno_notify_acceptable_iterate_callba } bool uno_set_logger_stream_callback(uno_logger_stream_callback logger_stream_callback, void* user_data) { - delete c_ostream; - c_ostream = new COStream(logger_stream_callback, user_data); - Logger::set_stream(*c_ostream); + delete c_stream_callback; + delete ostream; + c_stream_callback = new CStreamCallback(logger_stream_callback, user_data); + ostream = new UserOStream(c_stream_callback); + Logger::set_stream(*ostream); return true; } bool uno_reset_logger_stream() { - delete c_ostream; - c_ostream = nullptr; + delete ostream; + delete c_stream_callback; + c_stream_callback = nullptr; + ostream = nullptr; Logger::set_stream(std::cout); return true; } diff --git a/interfaces/Matlab/README.md b/interfaces/Matlab/README.md new file mode 100644 index 000000000..80b0cc5b1 --- /dev/null +++ b/interfaces/Matlab/README.md @@ -0,0 +1,272 @@ +# Uno's Matlab interface + +Uno's Matlab interface allows you to solve an optimization model described in Matlab using UNO, a modern and modular solver for nonlinearly constrained optimization. + +The package can be used in two ways: + +- high-level fmincon-style interface through the Matlab function `uno`; this is the simplest way to employ Uno from Matlab; +- low-level interface through the Matlab MEX functions `uno_optimize` and `uno_options`; this gives more control with advanced solver features. + +An example is available in the file [example_hs015.m](example/example_hs015.m). + +## Affiliation + +This Matlab interface is developed and maintained by [Stefano Lovato](https://github.com/stefphd) and [Charlie Vanaret](https://github.com/cvanaret). + +## Installation + +To install UNO's Matlab interface add the package directory to the Matlab path. The directory is either `/bin` or `/lib`, depending on your system. + +## High-level interface + +The high-level interface provides a fmincon-style syntax for defining the optimization problem and calling the UNO solver. It runs on top of the low-level interface, but note that not all UNO features are available this way. + +Use `help uno` to display the help text for `uno`. + +Define the objective function and its gradient: + +```matlab +function [fval,grad] = fun(x) + % ... +end +``` + +- the objective value `fval` must be a scalar numeric value; +- the objective gradient `grad` must be a numeric vector of `length(x)` elements; + +Define the nonlinear constraint function `c(x)<=0` and `ceq(x)=0` and its gradient: + +```matlab +function [c,ceq,gradc,gradceq] = nlcon(x) + % ... +end +``` + +- the inequality `c` and equality `ceq` constraints must be a numeric vector (possibly empty for no constraint); +- the inequality `gradc` and equality `gradceq` constraint gradients must be numeric matrices with sizes `length(x)`-by-`length(c)` and `length(x)`-by-`length(ceq)`, respectively. Each column is the gradient of a constraint (possibly empty for no constraint); + +Define the linear constraints `A*x<=b` and `Aeq*x=beq`: + +```matlab +A = [ ... ]; +b = [ ... ]; +Aeq = [ ... ]; +beq = [ ... ]; +``` + +- `b` and `beq` must be numeric vector (possibly empty for no constraint); +- `A` and `Aeq` must be numeric matrices with sizes `length(b)`-by-`length(x)` and `length(beq)`-by-`length(x)`, respectively (possibly empty for no constraint); + +Define the variable lower and upper bounds: + +```matlab +lb = [ ... ]; +ub = [ ... ]; +``` + +Define the initial guess: + +```matlab +x0 = [ ... ]; +``` + +Optionally define the options: + +```matlab +options = struct(); +% ... +``` + +The default options can be obtained: + +```matlab +options = uno('defaults') +% or options = uno('defaults',preset) to get the options for a given preset +``` + +The low-level interface MEX function `uno_options` can be also used to obtain the default options. + +Additionaly, define the Lagrangian Hessian function in the options as (similarly to Matlab fmincon): + +```matlab +options.HessianFnc = @(x,rho,lambda) hessian_fnc(x,rho,lambda); + +function H = hessian_fnc(x,rho,lambda) + % ... +end +``` + +- `rho` is the objective multiplier; +- `lambda.ineqnonlin` and `lambda.eqnonlin` are the Lagrange multipliers; +- `H` must be a symmetric numeric matrix with size `length(x)`-by-`length(x)` + +Finally, call `uno` function using a fmincon-style syntax: + +``` +[x,fval,exitflag,output,lambda,grad,hessian] = uno(@fun,x0,A,b,Aeq,beq,lb,ub,@nlcon,options); +``` + +The major differences between `uno` and `fmincon` syntax are that: + +- objective and nonlinear constraint gradients must be always given (no finite difference is available); to employ finite difference, users must implement them explicitly; +- the objective multiplier is also required in the Lagrangian Hessian function; if not provided, LBFGS is employed. + +## Low-level interface + +The low-level interface gives you full access to all UNO options and functionality. Use it if you need more control with advanced solver features. + +Use `help uno_options` and `help uno_optimize` to display the help text for `uno_options` and `uno_optimize`. + +### Building an optimization model + +Create optimization model as a Matlab struct: + +```matlab +% Problem type: 'L' = Linear, 'Q' = Quadratic, 'N' = Nonlinear +model.problem_type = ...; +% Vector indexing: 1 = Matlab-style, 0 = C-style +model.base_indexing = ...; +% Number of variables +model.number_variables = ...; +% Variable bounds +model.variables_lower_bounds = [ ... ]; +model.variables_upper_bounds = [ ... ]; +``` + +The following fields has also to be set: + +- objective and its full gradient; +```matlab +% Optimization sense: 1 = Minimize, -1 = Maximize +model.optimization_sense = 1; +% Objective function handle: objective = objective_function(x) +model.objective_function = @objective_function; +% Objective gradient handle: gradient = objective_gradient(x) +model.objective_gradient = @objective_gradient; +``` + +- constraint function and bounds and its sparse Jacobian matrix; +```matlab +% Number of constraints +model.number_constraints = ...; +% Constraint bounds +model.constraints_lower_bounds = [ ... ]; +model.constraints_upper_bounds = [ ... ]; +% Constraint function handle: constraints = constrain_function(x) +model.constraint_function = @constraint_function; +% Constraint jacobian handle: jacobian_values = constraint_jacobian(x) +model.constraint_jacobian = @constraint_jacobian; +% Constraint sparsity pattern (base_indexing-based) +model.number_jacobian_nonzeros = ...; +model.jacobian_row_indices = [ ... ]; +model.jacobian_column_indices = [ ... ]; +``` + +- sparse Lagrangian Hessian matrix; +```matlab +% Lagrangian sign convention: 1 = rho*f(x) + y^T c(x), -1 = rho*f(x) - y^T c(x) +model.lagrangian_sign_convention = ...; +% Hessian triangular part: 'L' = lower, 'U' = upper +model.hessian_triangular_part = 'L'; +% Lagrangian Hessian handle: hessian_values = lagrangian_hessian(x,rho,y) +model.lagrangian_hessian = @lagrangian_hessian; +% Lagrangian Hessian sparsity pattern (base_indexing-based) +model.number_hessian_nonzeros = ...; +model.hessian_row_indices = [ ... ]; +model.hessian_column_indices = [ ... ]; +``` + +- a Jacobian operator (performs Jacobian-vector products); +```matlab +% Jacobian operator handle: result = jacobian_operator(x,v) +model.jacobian_operator = @jacobian_operator; +``` + +- a Jacobian-transposed operator (performs Jacobian-transposed-vector products); +```matlab +% Jacobian tranposed operator handle: result = jacobian_transposed_operator(x,v) +model.jacobian_transposed_operator = @jacobian_transposed_operator; +``` + +- a Hessian operator (performs Hessian-vector products); +```matlab +% Hessian operator handle: result = lagrangian_hessian_operator(x,rho,y,v) +model.lagrangian_hessian_operator = @lagrangian_hessian_operator; +``` + +- an initial primal and dual point; +```matlab +% Initial primal point +model.initial_primal_iterate = [ ... ]; +% Dual primal point +model.initial_dual_iterate = [ ... ]; +``` + +User data can be employed in the optimization model functions using function handle, for example +```matlab +model.objective_function = @(x) objective_function(x, user_data); +``` + +### Defining UNO solver options + +Optionally obtain the default options for the UNO solver as a Matlab struct: +```matlab +options = uno_options(); +% or options = uno_options(preset) to get the options for a given preset +``` + +Default preset is used if empty or no arguments are given + +### Defining UNO solver callbacks + +Optionally define UNO solver callbacks in a Matlab struct: + +- Logger stream callback (called when printing string to the stream) +```matlab +% Logger stream handle: logger_stream_callback(str) +callbacks.logger_stream_callback = @logger_stream_callback; +``` + +- Notify callback (called when acceptable iterate found) +```matlab +% Notify acceptable iterate handle: notify_acceptable_iterate_callback(x, yl, yb, y, rho, feas, stat, compl) +callbacks.notify_acceptable_iterate_callback = @notify_acceptable_iterate_callback; +``` + +- User termination callback (called at the major iteration) +```matlab +% User termination handle: terminate = user_termination_callback(x, yl, yb, y, rho, feas, stat, compl) +callbacks.user_termination_callback = @user_termination_callback; +``` + +### Solving the model + +The model can then be solved by Uno, with optional input arguments `options` and `callbacks`: +```matlab +% Call UNO optimizer +% result = uno_optimize(model); +% result = uno_optimize(model, options); +result = uno_optimize(model, options, callbacks); +``` + +### Inspecting the result + +To inspect the result of the optimization, read the fields of the `result` struct: +- the optimization status (`UNO_SUCCESS = 0`, `UNO_ITERATION_LIMIT = 1`, `UNO_TIME_LIMIT = 2`, `UNO_EVALUATION_ERROR = 3`, `UNO_ALGORITHMIC_ERROR = 4`, `UNO_USER_TERMINATION = 5`): `result.optimization_status` +- the solution status (`UNO_NOT_OPTIMAL = 0`, `UNO_FEASIBLE_KKT_POINT = 1`, `UNO_FEASIBLE_FJ_POINT = 2`, `UNO_INFEASIBLE_STATIONARY_POINT = 3`, `UNO_FEASIBLE_SMALL_STEP = 4`, `UNO_INFEASIBLE_SMALL_STEP = 5`, `UNO_UNBOUNDED = 6`): `result.solution_status` +- the objective value of the solution: `result.solution_objective` +- the primal solution: `result.primal_solution` +- the dual solution associated with the general constraints: `result.constraint_dual_solution` +- the dual solution associated with the lower bounds: `result.lower_bound_dual_solution` +- the dual solution associated with the upper bounds: `result.upper_bound_dual_solution` +- the primal feasibility residual at the solution: `result.solution_primal_feasibility` +- the stationarity residual at the solution: `result.solution_stationarity` +- the complementarity residual at the solution: `result.solution_complementarity` +- the number of (outer) iterations: `result.number_iterations` +- the CPU time (in seconds): `result.cpu_time` +- the number of objective evaluations: `result.number_objective_evaluations` +- the number of constraint evaluations: `result.number_constraint_evaluations` +- the number of objective gradient evaluations: `result.number_objective_gradient_evaluations` +- the number of Jacobian evaluations: `result.number_jacobian_evaluations` +- the number of Hessian evaluations: `result.number_hessian_evaluations` +- the number of subproblems solved: `result.number_subproblems_solved` \ No newline at end of file diff --git a/interfaces/Matlab/cpp_classes/ErrorString.cpp b/interfaces/Matlab/cpp_classes/ErrorString.cpp new file mode 100644 index 000000000..ca132d7d5 --- /dev/null +++ b/interfaces/Matlab/cpp_classes/ErrorString.cpp @@ -0,0 +1,51 @@ +// Copyright (c) 2025 Stefano Lovato and Charlie Vanaret +// Licensed under the MIT license. See LICENSE file in the project directory for details. + +#include "ErrorString.hpp" +#include + +namespace uno { + + const std::map ErrorString::error_strings = { + {ErrorType::NARGIN_NOTENOUGH, "Invalid argument list. Function requires %d more input(s)."}, + {ErrorType::NARGIN_TOOMANY, "Too many input arguments."}, + {ErrorType::INPUT_INVALID, "Invalid input argument at position %d."}, + {ErrorType::INPUT_STRING, "Invalid argument at position %d. Value must be char or string."}, + {ErrorType::INPUT_STRUCT, "Invalid argument at position %d. Value must be struct."}, + {ErrorType::INPUT_SCALAR, "Invalid input argument at position %d. Value must be scalar."}, + {ErrorType::OUTPUT_SCALAR, "Invalid output argument at position %d. Value must be a numeric scalar."}, + {ErrorType::OUTPUT_VECTOR, "Invalid output argument at position %d. Value must be a numeric vector of %d element(s)."}, + {ErrorType::OUTPUT_BOOL_SCALAR, "Invalid output argument at position %d. Value must be a logical scalar."}, + {ErrorType::INVALID_HANDLE, "Value must be a function handle."}, + {ErrorType::EVALUATION, "Error in function evaluation."}, + {ErrorType::OPTION_TYPE, "Invalid type of option '%s'."}, + {ErrorType::MISSING_FIELD, "Missing field '%s'."}, + {ErrorType::FIELD_CHAR, "Type of field '%s' must be char."}, + {ErrorType::FIELD_STRING, "Type of field '%s' must be char array or string."}, + {ErrorType::FIELD_POSITIVE_INTEGER, "Field '%s' must be a positive integer value."}, + {ErrorType::FIELD_UNITARY, "Field '%s' must be either +1 or -1."}, + {ErrorType::FIELD_VECTOR, "Field '%s' must be a vector of %d element(s)."}, + {ErrorType::FIELD_VECTOR_INT, "Field '%s' must be a vector of %d integer element(s)."}, + {ErrorType::UNO, "Error in UNO. %s"} + }; + + std::string ErrorString::format_error(const ErrorType error_type, ...) { + const auto it = error_strings.find(error_type); + if (it == error_strings.end()) { + return "Unknown error."; + } + // format the string (avoid using std::format which is C++20) + const std::string& error_fmt = it->second; + va_list args1; + va_start(args1, error_type); + va_list args2; + va_copy(args2, args1); + const int size = std::vsnprintf(nullptr, 0, error_fmt.c_str(), args2); + va_end(args2); + std::string error_msg(size, '\0'); + std::vsnprintf(error_msg.data(), error_msg.size() + 1, error_fmt.c_str(), args1); + va_end(args1); + return error_msg; + } + +} // namespace \ No newline at end of file diff --git a/interfaces/Matlab/cpp_classes/ErrorString.hpp b/interfaces/Matlab/cpp_classes/ErrorString.hpp new file mode 100644 index 000000000..e747c3133 --- /dev/null +++ b/interfaces/Matlab/cpp_classes/ErrorString.hpp @@ -0,0 +1,45 @@ +// Copyright (c) 2025 Stefano Lovato and Charlie Vanaret +// Licensed under the MIT license. See LICENSE file in the project directory for details. + +#ifndef UNO_ERRORSTRING_H +#define UNO_ERRORSTRING_H + +#include +#include + +namespace uno { + + enum class ErrorType { + NARGIN_NOTENOUGH, + NARGIN_TOOMANY, + INPUT_INVALID, + INPUT_STRING, + INPUT_STRUCT, + INPUT_SCALAR, + OUTPUT_SCALAR, + OUTPUT_VECTOR, + OUTPUT_BOOL_SCALAR, + INVALID_HANDLE, + EVALUATION, + OPTION_TYPE, + MISSING_FIELD, + FIELD_CHAR, + FIELD_STRING, + FIELD_POSITIVE_INTEGER, + FIELD_UNITARY, + FIELD_VECTOR, + FIELD_VECTOR_INT, + UNO + }; + + // map the error type to error strings + class ErrorString { + public: + static std::string format_error(ErrorType error_type, ...); + private: + static const std::map error_strings; + }; + +}; // namespace + +#endif // UNO_ERRORSTRING_H \ No newline at end of file diff --git a/interfaces/Matlab/cpp_classes/MatlabModel.cpp b/interfaces/Matlab/cpp_classes/MatlabModel.cpp new file mode 100644 index 000000000..917c434b7 --- /dev/null +++ b/interfaces/Matlab/cpp_classes/MatlabModel.cpp @@ -0,0 +1,400 @@ +// Copyright (c) 2025 Stefano Lovato and Charlie Vanaret +// Licensed under the MIT license. See LICENSE file in the project directory for details. + +#include +#include "MatlabModel.hpp" +#include "../unomex/unomex_function.hpp" +#include "../unomex/unomex_conversion.hpp" +#include "../unomex/unomex_validation.hpp" +#include "../unomex/unomex_utils.hpp" +#include "mex.h" + +namespace uno { + + MatlabModel::MatlabModel(const MatlabUserModel& user_model): + Model("Matlab model", static_cast(user_model.number_variables), static_cast(user_model.number_constraints), + static_cast(user_model.optimization_sense), user_model.lagrangian_sign_convention, user_model.base_indexing), + user_model(user_model), + nonlinear_constraints(this->number_constraints), + equality_constraints_collection(this->equality_constraints), + inequality_constraints_collection(this->inequality_constraints) { + this->find_fixed_variables(this->fixed_variables); + this->partition_constraints(this->equality_constraints, this->inequality_constraints); + } + + ProblemType MatlabModel::get_problem_type() const { + return this->user_model.problem_type; + } + + bool MatlabModel::has_jacobian_operator() const { + return this->user_model.jacobian_operator != nullptr; + } + + bool MatlabModel::has_jacobian_transposed_operator() const { + return this->user_model.jacobian_transposed_operator != nullptr; + } + + bool MatlabModel::has_hessian_operator() const { + return this->user_model.lagrangian_hessian_operator != nullptr; + } + + bool MatlabModel::has_hessian_matrix() const { + return this->user_model.lagrangian_hessian != nullptr; + } + + double MatlabModel::evaluate_objective(const Vector& x) const { + double objective_value{0.}; + // objective_value = objective_function(x); + if (this->user_model.objective_function) { + const std::vector inputs({vector_to_mxArray(x, this->number_variables)}); + std::vector outputs(1); + MxArrayVectorGuard input_guard(inputs); + MxArrayVectorGuard output_guard(outputs); + try { + call_matlab_function(user_model.objective_function, inputs, outputs); + } catch (const MatlabFunctionError& err) { + mexWarnMsgIdAndTxt("uno:warning", "%s\n", err.what()); + throw FunctionEvaluationError(); + } + if (std::string errmsg; !validate_double_scalar_output(outputs[0], errmsg)) { + mexWarnMsgIdAndTxt("uno:warning", "Error in objective function.\n%s\n", errmsg.c_str()); + throw FunctionEvaluationError(); + } + objective_value = this->optimization_sense * mxArray_to_scalar(outputs[0]); + ++this->number_model_evaluations.objective; + } + return objective_value; + } + + void MatlabModel::evaluate_constraints(const Vector& x, Vector& constraints) const { + // constraint = constraint_function(x); + if (this->user_model.constraint_functions) { + std::vector inputs({vector_to_mxArray(x, this->number_variables)}); + std::vector outputs(1); + MxArrayVectorGuard input_guard(inputs); + MxArrayVectorGuard output_guard(outputs); + try { + call_matlab_function(user_model.constraint_functions, inputs, outputs); + } catch (const MatlabFunctionError& err) { + mexWarnMsgIdAndTxt("uno:warning", "%s", err.what()); + throw FunctionEvaluationError(); + } + if (std::string errmsg; !validate_double_vector_output(outputs[0], this->number_constraints, errmsg)) { + mexWarnMsgIdAndTxt("uno:warning", "Error in constraint function.\n%s\n", errmsg.c_str()); + throw FunctionEvaluationError(); + } + mxArray_to_vector(outputs[0], constraints); + ++this->number_model_evaluations.constraints; + } + } + + void MatlabModel::evaluate_objective_gradient(const Vector& x, Vector& gradient) const { + // gradient = objective_gradient(x); + if (this->user_model.objective_gradient) { + std::vector inputs({vector_to_mxArray(x, this->number_variables)}); + std::vector outputs(1); + MxArrayVectorGuard input_guard(inputs); + MxArrayVectorGuard output_guard(outputs); + try { + call_matlab_function(user_model.objective_gradient, inputs, outputs); + } catch (const MatlabFunctionError& err) { + mexWarnMsgIdAndTxt("uno:warning", "%s", err.what()); + throw GradientEvaluationError(); + } + if (std::string errmsg; !validate_double_vector_output(outputs[0], this->number_variables, errmsg)) { + mexWarnMsgIdAndTxt("uno:warning", "Error in objective gradient.\n%s\n", errmsg.c_str()); + throw FunctionEvaluationError(); + } + mxArray_to_vector(outputs[0], gradient); + for (size_t variable_index: Range(this->number_variables)) { + gradient[variable_index] *= this->optimization_sense; + } + ++this->number_model_evaluations.objective_gradient; + } + } + + void MatlabModel::compute_jacobian_sparsity(uno_int* row_indices, uno_int* column_indices, uno_int row_offset, uno_int column_offset, + uno_int solver_indexing, MatrixOrder /*matrix_order*/) const { + // copy the indices of the user sparsity patterns to the Uno vectors + for (size_t nonzero_index: Range(static_cast(this->user_model.number_jacobian_nonzeros))) { + row_indices[nonzero_index] = this->user_model.jacobian_row_indices[nonzero_index] + row_offset; + column_indices[nonzero_index] = this->user_model.jacobian_column_indices[nonzero_index] + column_offset; + } + // TODO matrix_order + + // handle the solver indexing + if (this->user_model.base_indexing != solver_indexing) { + const int indexing_difference = solver_indexing - this->user_model.base_indexing; + for (size_t nonzero_index: Range(static_cast(this->user_model.number_jacobian_nonzeros))) { + row_indices[nonzero_index] += indexing_difference; + column_indices[nonzero_index] += indexing_difference; + } + } + } + + void MatlabModel::compute_hessian_sparsity(uno_int* row_indices, uno_int* column_indices, uno_int solver_indexing) const { + // copy the indices of the user sparsity patterns to the Uno vectors + const size_t number_hessian_nonzeros = this->number_hessian_nonzeros(); + for (size_t nonzero_index: Range(number_hessian_nonzeros)) { + row_indices[nonzero_index] = this->user_model.hessian_row_indices[nonzero_index]; + column_indices[nonzero_index] = this->user_model.hessian_column_indices[nonzero_index]; + } + + // handle the solver indexing + if (this->user_model.base_indexing != solver_indexing) { + const int indexing_difference = solver_indexing - this->user_model.base_indexing; + for (size_t nonzero_index: Range(number_hessian_nonzeros)) { + row_indices[nonzero_index] += indexing_difference; + column_indices[nonzero_index] += indexing_difference; + } + } + } + + void MatlabModel::evaluate_jacobian(const Vector& x, double* jacobian_values) const { + // jacobian_values = jacobian(x); + if (this->user_model.jacobian) { + std::vector inputs({vector_to_mxArray(x, this->number_variables)}); + std::vector outputs(1); + MxArrayVectorGuard input_guard(inputs); + MxArrayVectorGuard output_guard(outputs); + try { + call_matlab_function(user_model.jacobian, inputs, outputs); + } catch (const MatlabFunctionError& err) { + mexWarnMsgIdAndTxt("uno:warning", "%s", err.what()); + throw GradientEvaluationError(); + } + if (std::string errmsg; !validate_double_vector_output(outputs[0], this->number_jacobian_nonzeros(), errmsg)) { + mexWarnMsgIdAndTxt("uno:warning", "Error in Jacobian.\n%s\n", errmsg.c_str()); + throw FunctionEvaluationError(); + } + mxArray_to_pointer(outputs[0], jacobian_values); + ++this->number_model_evaluations.jacobian; + } + } + + void MatlabModel::evaluate_lagrangian_hessian(const Vector& x, double objective_multiplier, const Vector& multipliers, + double* hessian_values) const { + // hessian_values = lagrangian_hessian(x, rho, y); + if (this->user_model.lagrangian_hessian) { + objective_multiplier *= this->optimization_sense; + // if the model has a different sign convention for the Lagrangian than Uno, flip the signs of the multipliers + if (this->user_model.lagrangian_sign_convention == UNO_MULTIPLIER_POSITIVE) { + const_cast&>(multipliers).scale(-1.); + } + // eval matlab function + std::vector inputs({vector_to_mxArray(x, this->number_variables), scalar_to_mxArray(objective_multiplier),vector_to_mxArray(multipliers)}); + std::vector outputs(1); + MxArrayVectorGuard input_guard(inputs); + MxArrayVectorGuard output_guard(outputs); + // flip the signs of the multipliers back + if (this->user_model.lagrangian_sign_convention == UNO_MULTIPLIER_POSITIVE) { + const_cast&>(multipliers).scale(-1.); + } + try { + call_matlab_function(user_model.lagrangian_hessian, inputs, outputs); + } catch (const MatlabFunctionError& err) { + mexWarnMsgIdAndTxt("uno:warning", "%s", err.what()); + throw HessianEvaluationError(); + } + if (std::string errmsg; !validate_double_vector_output(outputs[0], this->number_hessian_nonzeros(), errmsg)) { + mexWarnMsgIdAndTxt("uno:warning", "Error in Lagrangian Hessian.\n%s\n", errmsg.c_str()); + throw FunctionEvaluationError(); + } + mxArray_to_pointer(outputs[0], hessian_values); + ++this->number_model_evaluations.hessian; + } + else { + throw std::runtime_error("evaluate_lagrangian_hessian not implemented"); + } + } + + void MatlabModel::compute_jacobian_vector_product(const double* x, const double* vector, double* result) const { + // result = jacobian_operator(x, v); + if (this->user_model.jacobian_operator) { + mxArray* x_arr = pointer_to_mxArray(x, this->number_variables); + mxArray* vector_arr = pointer_to_mxArray(vector, this->number_variables); + std::vector inputs({x_arr, vector_arr}); + std::vector outputs(1); + MxArrayVectorGuard input_guard(inputs); + MxArrayVectorGuard output_guard(outputs); + try { + call_matlab_function(user_model.jacobian_operator, inputs, outputs); + } catch (const MatlabFunctionError& err) { + mexWarnMsgIdAndTxt("uno:warning", "%s", err.what()); + throw GradientEvaluationError(); + } + if (std::string errmsg; !validate_double_vector_output(outputs[0], this->number_constraints, errmsg)) { + mexWarnMsgIdAndTxt("uno:warning", "Error in Jacobian operator.\n%s\n", errmsg.c_str()); + throw FunctionEvaluationError(); + } + mxArray_to_pointer(outputs[0], result); + } + else { + throw std::runtime_error("compute_jacobian_vector_product not implemented"); + } + } + + void MatlabModel::compute_jacobian_transposed_vector_product(const double* x, const double* vector, double* result) const { + // result = jacobian_transposed_operator(x, v); + if ((this->user_model.jacobian_operator != nullptr) && + !mxIsEmpty(this->user_model.jacobian_operator)) { + mxArray* x_arr = pointer_to_mxArray(x, this->number_variables); + mxArray* vector_arr = pointer_to_mxArray(vector, this->number_constraints); + std::vector inputs({x_arr, vector_arr}); + std::vector outputs(1); + MxArrayVectorGuard input_guard(inputs); + MxArrayVectorGuard output_guard(outputs); + try { + call_matlab_function(user_model.jacobian_transposed_operator, inputs, outputs); + } catch (const MatlabFunctionError& err) { + mexWarnMsgIdAndTxt("uno:warning", "%s", err.what()); + throw GradientEvaluationError(); + } + if (std::string errmsg; !validate_double_vector_output(outputs[0], this->number_variables, errmsg)) { + mexWarnMsgIdAndTxt("uno:warning", "Error in Jacobian transposed operator.\n%s\n", errmsg.c_str()); + throw FunctionEvaluationError(); + } + mxArray_to_pointer(outputs[0], result); + } + else { + throw std::runtime_error("compute_jacobian_transposed_vector_product not implemented"); + } + } + + void MatlabModel::compute_hessian_vector_product(const double* x, const double* vector, double objective_multiplier, const Vector& multipliers, + double* result) const { + // result = lagrangian_hessian_operator(x, rho, y, v); + if (this->user_model.lagrangian_hessian_operator) { + objective_multiplier *= this->optimization_sense; + // if the model has a different sign convention for the Lagrangian than Uno, flip the signs of the multipliers + if (this->user_model.lagrangian_sign_convention == UNO_MULTIPLIER_POSITIVE) { + const_cast&>(multipliers).scale(-1.); + } + mxArray* x_arr = pointer_to_mxArray(x, this->number_variables); + mxArray* vector_arr = pointer_to_mxArray(vector, this->number_variables); + std::vector inputs({x_arr, scalar_to_mxArray(objective_multiplier), vector_to_mxArray(multipliers), vector_arr}); + std::vector outputs(1); + MxArrayVectorGuard input_guard(inputs); + MxArrayVectorGuard output_guard(outputs); + // flip the signs of the multipliers back + if (this->user_model.lagrangian_sign_convention == UNO_MULTIPLIER_POSITIVE) { + const_cast&>(multipliers).scale(-1.); + } + try { + call_matlab_function(user_model.lagrangian_hessian_operator, inputs, outputs); + } catch (const MatlabFunctionError& err) { + mexWarnMsgIdAndTxt("uno:warning", "%s", err.what()); + throw HessianEvaluationError(); + } + if (std::string errmsg; !validate_double_vector_output(outputs[0], this->number_variables, errmsg)) { + mexWarnMsgIdAndTxt("uno:warning", "Error in Lagrangian Hessian operator.\n%s\n", errmsg.c_str()); + throw FunctionEvaluationError(); + } + mxArray_to_pointer(outputs[0], result); + } + else { + throw std::runtime_error("compute_hessian_vector_product not implemented"); + } + } + + const std::vector& MatlabModel::get_variables_lower_bounds() const { + return this->user_model.variables_lower_bounds; + } + + const std::vector& MatlabModel::get_variables_upper_bounds() const { + return this->user_model.variables_upper_bounds; + } + + const SparseVector& MatlabModel::get_slacks() const { + return this->slacks; + } + + const Vector& MatlabModel::get_fixed_variables() const { + return this->fixed_variables; + } + + const std::vector& MatlabModel::get_constraints_lower_bounds() const { + return this->user_model.constraints_lower_bounds; + } + + const std::vector& MatlabModel::get_constraints_upper_bounds() const { + return this->user_model.constraints_upper_bounds; + } + + const Collection& MatlabModel::get_equality_constraints() const { + return this->equality_constraints_collection; + } + + const Collection& MatlabModel::get_inequality_constraints() const { + return this->inequality_constraints_collection; + } + + const Collection& MatlabModel::get_linear_constraints() const { + return this->linear_constraints; + } + + const Collection& MatlabModel::get_nonlinear_constraints() const { + return this->nonlinear_constraints; + } + + void MatlabModel::initial_primal_point(Vector& x) const { + // copy the initial primal point + for (const size_t variable_index: Range(static_cast(this->user_model.number_variables))) { + x[variable_index] = this->user_model.initial_primal_iterate[variable_index]; + } + } + + void MatlabModel::initial_dual_point(Vector& multipliers) const { + // copy the initial dual point + for (const size_t constraint_index: Range(static_cast(this->user_model.number_constraints))) { + multipliers[constraint_index] = this->user_model.initial_dual_iterate[constraint_index]; + } + if (this->user_model.lagrangian_sign_convention == UNO_MULTIPLIER_POSITIVE) { + multipliers.scale(-1.); + } + } + + void MatlabModel::postprocess_solution(Iterate& /* iterate */) const { + // do nothing + } + + size_t MatlabModel::number_jacobian_nonzeros() const { + return static_cast(this->user_model.number_jacobian_nonzeros); + } + + size_t MatlabModel::number_hessian_nonzeros() const { + if (this->user_model.number_hessian_nonzeros.has_value()) { + return static_cast(*this->user_model.number_hessian_nonzeros); + } + else { + throw std::runtime_error("The number of Hessian nonzeros is not available in MatlabModel"); + } + } + + size_t MatlabModel::number_model_objective_evaluations() const { + return this->number_model_evaluations.objective; + } + + size_t MatlabModel::number_model_constraints_evaluations() const { + return this->number_model_evaluations.constraints; + } + + size_t MatlabModel::number_model_objective_gradient_evaluations() const { + return this->number_model_evaluations.objective_gradient; + } + + size_t MatlabModel::number_model_jacobian_evaluations() const { + return this->number_model_evaluations.jacobian; + } + + size_t MatlabModel::number_model_hessian_evaluations() const { + return this->number_model_evaluations.hessian; + } + + void MatlabModel::reset_number_evaluations() const { + this->number_model_evaluations.reset(); + } + +} // namespace + + diff --git a/interfaces/Matlab/cpp_classes/MatlabModel.hpp b/interfaces/Matlab/cpp_classes/MatlabModel.hpp new file mode 100644 index 000000000..3a0828860 --- /dev/null +++ b/interfaces/Matlab/cpp_classes/MatlabModel.hpp @@ -0,0 +1,103 @@ +// Copyright (c) 2025 Stefano Lovato and Charlie Vanaret +// Licensed under the MIT license. See LICENSE file in the project directory for details. + +#ifndef UNO_MATLABMODEL_H +#define UNO_MATLABMODEL_H + +#include "Uno.hpp" +#include "linear_algebra/SparseVector.hpp" +#include "linear_algebra/Vector.hpp" +#include "model/Model.hpp" +#include "symbolic/IntegerRange.hpp" +#include "symbolic/CollectionAdapter.hpp" +#include "tools/NumberModelEvaluations.hpp" +#include "../../UserModel.hpp" +#include "../unomex/unomex_function.hpp" +#include "mex.h" + +namespace uno { + + using MatlabUserModel = UserModel, void*>; // user data are handled in matlab wrapper + + // MatlabModel contains an instance of UserModel and complies with the Model interface + class MatlabModel : public Model { + public: + MatlabModel(const MatlabUserModel& user_model); + ~MatlabModel() override = default; + + [[nodiscard]] ProblemType get_problem_type() const override; + + // availability of linear operators + [[nodiscard]] bool has_jacobian_operator() const override; + [[nodiscard]] bool has_jacobian_transposed_operator() const override; + [[nodiscard]] bool has_hessian_operator() const override; + [[nodiscard]] bool has_hessian_matrix() const override; + + // function evaluations + [[nodiscard]] double evaluate_objective(const Vector& x) const override; + void evaluate_constraints(const Vector& x, Vector& constraints) const override; + + // dense objective gradient + void evaluate_objective_gradient(const Vector& x, Vector& gradient) const override; + + // sparsity patterns of Jacobian and Hessian + void compute_jacobian_sparsity(uno_int* row_indices, uno_int* column_indices, uno_int row_offset, uno_int column_offset, + uno_int solver_indexing, MatrixOrder /*matrix_order*/) const override; + void compute_hessian_sparsity(uno_int* row_indices, uno_int* column_indices, uno_int solver_indexing) const override; + + // numerical evaluations of Jacobian and Hessian + void evaluate_jacobian(const Vector& x, double* jacobian_values) const override; + void evaluate_lagrangian_hessian(const Vector& x, double objective_multiplier, const Vector& multipliers, + double* hessian_values) const override; + + // linear operators for Jacobian-, Jacobian^T-, and Hessian-vector products + void compute_jacobian_vector_product(const double* x, const double* vector, double* result) const override; + void compute_jacobian_transposed_vector_product(const double* x, const double* vector, double* result) const override; + void compute_hessian_vector_product(const double* x, const double* vector, double objective_multiplier, const Vector& multipliers, + double* result) const override; + + // purely functions + [[nodiscard]] const std::vector& get_variables_lower_bounds() const override; + [[nodiscard]] const std::vector& get_variables_upper_bounds() const override; + [[nodiscard]] const SparseVector& get_slacks() const override; + [[nodiscard]] const Vector& get_fixed_variables() const override; + + [[nodiscard]] const std::vector& get_constraints_lower_bounds() const override; + [[nodiscard]] const std::vector& get_constraints_upper_bounds() const override; + + [[nodiscard]] const Collection& get_equality_constraints() const override; + [[nodiscard]] const Collection& get_inequality_constraints() const override; + [[nodiscard]] const Collection& get_linear_constraints() const override; + [[nodiscard]] const Collection& get_nonlinear_constraints() const override; + + void initial_primal_point(Vector& x) const override; + void initial_dual_point(Vector& multipliers) const override; + void postprocess_solution(Iterate& iterate) const override; + + [[nodiscard]] size_t number_jacobian_nonzeros() const override; + [[nodiscard]] size_t number_hessian_nonzeros() const override; + + [[nodiscard]] size_t number_model_objective_evaluations() const override; + [[nodiscard]] size_t number_model_constraints_evaluations() const override ; + [[nodiscard]] size_t number_model_objective_gradient_evaluations() const override; + [[nodiscard]] size_t number_model_jacobian_evaluations() const override; + [[nodiscard]] size_t number_model_hessian_evaluations() const override; + void reset_number_evaluations() const override; + + protected: + const MatlabUserModel& user_model; + mutable NumberModelEvaluations number_model_evaluations; + const SparseVector slacks{}; + Vector fixed_variables{}; + const IntegerRange linear_constraints{0}; + const IntegerRange nonlinear_constraints; + std::vector equality_constraints; + CollectionAdapter> equality_constraints_collection; + std::vector inequality_constraints; + CollectionAdapter> inequality_constraints_collection; + }; + +}; // namespace + +#endif // UNO_MATLABMODEL_H \ No newline at end of file diff --git a/interfaces/Matlab/cpp_classes/MatlabStreamCallback.cpp b/interfaces/Matlab/cpp_classes/MatlabStreamCallback.cpp new file mode 100644 index 000000000..b64b57d66 --- /dev/null +++ b/interfaces/Matlab/cpp_classes/MatlabStreamCallback.cpp @@ -0,0 +1,47 @@ +// Copyright (c) 2025 Stefano Lovato and Charlie Vanaret +// Licensed under the MIT license. See LICENSE file in the project directory for details. + +#include +#include "MatlabStreamCallback.hpp" +#include "../unomex/unomex_function.hpp" +#include "../unomex/unomex_conversion.hpp" +#include "../unomex/unomex_validation.hpp" +#include "../unomex/unomex_utils.hpp" +#include "mex.h" + +namespace uno { + + MatlabStreamCallback::MatlabStreamCallback(handle_t logger_stream_callback) : + UserStreamCallback(), logger_stream_callback(logger_stream_callback) { } + + MatlabStreamCallback::~MatlabStreamCallback() { + mexEvalString("pause(0);"); // force output to appear immediately + } + + int32_t MatlabStreamCallback::operator()(const char* buf, int32_t len) const { + if (this->logger_stream_callback) { + // logger_stream_callback(str) + const std::string str(buf, len); + const std::vector inputs = { string_to_mxArray(str) }; + std::vector outputs(0); + MxArrayVectorGuard input_guard(inputs); + // MxArrayVectorGuard output_guard(outputs); + try { + call_matlab_function(this->logger_stream_callback, inputs, outputs); + } + catch (const MatlabFunctionError& err) { + mexWarnMsgIdAndTxt("uno:warning", "%s", err.what()); + } + // no output to validate + mexEvalString("pause(0);"); // force output to appear immediately + return len; + } + else { + // call mexPrintf + mexPrintf("%.*s", static_cast(len), buf); + mexEvalString("pause(0);"); // force output to appear immediately + return len; + } + } + +} // namespace diff --git a/interfaces/Matlab/cpp_classes/MatlabStreamCallback.hpp b/interfaces/Matlab/cpp_classes/MatlabStreamCallback.hpp new file mode 100644 index 000000000..aaf4c5e07 --- /dev/null +++ b/interfaces/Matlab/cpp_classes/MatlabStreamCallback.hpp @@ -0,0 +1,25 @@ +// Copyright (c) 2025 Stefano Lovato and Charlie Vanaret +// Licensed under the MIT license. See LICENSE file in the project directory for details. + +#ifndef UNO_MATLABSTREAMCALLBACK_H +#define UNO_MATLABSTREAMCALLBACK_H + +#include "../../UserModel.hpp" +#include "../unomex/unomex_function.hpp" + +namespace uno { + + class MatlabStreamCallback : public UserStreamCallback { + public: + explicit MatlabStreamCallback(handle_t logger_stream_callback); + ~MatlabStreamCallback() override; + + int32_t operator()(const char* buf, int32_t len) const override; + + private: + handle_t logger_stream_callback; + }; + +}; // namespace + +#endif // UNO_MATLABSTREAMCALLBACK_H \ No newline at end of file diff --git a/interfaces/Matlab/cpp_classes/MatlabUserCallbacks.cpp b/interfaces/Matlab/cpp_classes/MatlabUserCallbacks.cpp new file mode 100644 index 000000000..0340ace8a --- /dev/null +++ b/interfaces/Matlab/cpp_classes/MatlabUserCallbacks.cpp @@ -0,0 +1,75 @@ +// Copyright (c) 2025 Stefano Lovato and Charlie Vanaret +// Licensed under the MIT license. See LICENSE file in the project directory for details. + +#include +#include "MatlabUserCallbacks.hpp" +#include "../unomex/unomex_function.hpp" +#include "../unomex/unomex_conversion.hpp" +#include "../unomex/unomex_validation.hpp" +#include "../unomex/unomex_utils.hpp" +#include "mex.h" + +#ifdef HAS_MXLIBUT +// These functions are to allow interrupt from MATLAB using CTRL+C (see +// https://undocumentedmatlab.com/articles/mex-ctrl-c-interrupt) +extern "C" bool utIsInterruptPending(); +extern "C" bool utSetInterruptPending(bool); +#else +// libut is not available, use dummy +constexpr bool utIsInterruptPending() { return false; } +constexpr bool utSetInterruptPending(bool) { return false; } +#endif + +namespace uno { + + MatlabUserCallbacks::MatlabUserCallbacks(handle_t notify_acceptable_iterate_callback, handle_t termination_callback) : UserCallbacks(), + notify_acceptable_iterate_callback(notify_acceptable_iterate_callback), + termination_callback(termination_callback) { } + + void MatlabUserCallbacks::notify_acceptable_iterate(const Vector& primals, const Multipliers& multipliers, double objective_multiplier, double primal_feasibility, double stationarity, double complementarity) { + // notify_acceptable_iterate_callback(x, yl, yb, y, rho, feas, stat, compl); + if (this->notify_acceptable_iterate_callback) { + std::vector inputs = {vector_to_mxArray(primals), vector_to_mxArray(multipliers.lower_bounds), vector_to_mxArray(multipliers.upper_bounds), vector_to_mxArray(multipliers.constraints), scalar_to_mxArray(objective_multiplier), scalar_to_mxArray(primal_feasibility), scalar_to_mxArray(stationarity), scalar_to_mxArray(complementarity)}; + std::vector outputs(0); + MxArrayVectorGuard input_guard(inputs); + // MxArrayVectorGuard output_guard(outputs); + try { + call_matlab_function(this->notify_acceptable_iterate_callback, inputs, outputs); + } + catch (const MatlabFunctionError& err) { + mexWarnMsgIdAndTxt("uno:warning", "%s", err.what()); + } + // no output to validate + } + } + + bool MatlabUserCallbacks::termination(const Vector& primals, const Multipliers& multipliers, const double objective_multiplier, const double primal_feasibility, const double stationarity, const double complementarity) { + // handle matlab CTRL+C event + if (utIsInterruptPending()) { + utSetInterruptPending(false); + return true; + } + // terminate = user_termination_callback(x, yl, yb, y, rho, feas, stat, compl); + if (this->termination_callback) { + const std::vector inputs = {vector_to_mxArray(primals), vector_to_mxArray(multipliers.lower_bounds), vector_to_mxArray(multipliers.upper_bounds), vector_to_mxArray(multipliers.constraints), scalar_to_mxArray(objective_multiplier), scalar_to_mxArray(primal_feasibility), scalar_to_mxArray(stationarity), scalar_to_mxArray(complementarity)}; + std::vector outputs(1); + MxArrayVectorGuard input_guard(inputs); + MxArrayVectorGuard output_guard(outputs); + try { + call_matlab_function(this->termination_callback, inputs, outputs); + } + catch (const MatlabFunctionError& err) { + mexWarnMsgIdAndTxt("uno:warning", "%s", err.what()); + return false; + } + if (std::string errmsg; !validate_bool_scalar_output(outputs[0], errmsg)) { + mexWarnMsgIdAndTxt("uno:warning", "Error in user termination callback.\n%s\n", errmsg.c_str()); + return false; + } + return mxArray_to_scalar(outputs[0]); + } else { + return false; // never terminate + } + } + +} // namespace \ No newline at end of file diff --git a/interfaces/Matlab/cpp_classes/MatlabUserCallbacks.hpp b/interfaces/Matlab/cpp_classes/MatlabUserCallbacks.hpp new file mode 100644 index 000000000..9b55bfbab --- /dev/null +++ b/interfaces/Matlab/cpp_classes/MatlabUserCallbacks.hpp @@ -0,0 +1,28 @@ +// Copyright (c) 2025 Stefano Lovato and Charlie Vanaret +// Licensed under the MIT license. See LICENSE file in the project directory for details. + +#ifndef UNO_MATLABUSERCALLBACK_H +#define UNO_MATLABUSERCALLBACK_H + +#include "tools/UserCallbacks.hpp" +#include "../unomex/unomex_function.hpp" + +namespace uno { + + // Matlab user callbacks + class MatlabUserCallbacks : public UserCallbacks { + public: + MatlabUserCallbacks(handle_t notify_acceptable_iterate_callback, handle_t user_termination_callback); + + void notify_acceptable_iterate(const Vector& primals, const Multipliers& multipliers, double objective_multiplier, double primal_feasibility, double stationarity, double complementarity) override; + + bool termination(const Vector& primals, const Multipliers& multipliers, double objective_multiplier, double primal_feasibility, double stationarity, double complementarity) override; + + private: + handle_t notify_acceptable_iterate_callback; + handle_t termination_callback; + }; + +}; // namespace + +#endif // UNO_MATLABUSERCALLBACK_H \ No newline at end of file diff --git a/interfaces/Matlab/cpp_classes/MxStruct.cpp b/interfaces/Matlab/cpp_classes/MxStruct.cpp new file mode 100644 index 000000000..c752ebb8d --- /dev/null +++ b/interfaces/Matlab/cpp_classes/MxStruct.cpp @@ -0,0 +1,21 @@ +// Copyright (c) 2025 Stefano Lovato and Charlie Vanaret +// Licensed under the MIT license. See LICENSE file in the project directory for details. + +#include "MxStruct.hpp" + +namespace uno { + + mxArray* MxStruct::operator[](const std::string& key) const { + const auto it = this->data.find(key); + return (it != this->data.end()) ? it->second : nullptr; + } + + void MxStruct::insert(const std::string& key, mxArray* value) { + this->data[key] = value; + } + + bool MxStruct::contains(const std::string& key) const { + return this->data.find(key) != this->data.end(); + } + +} // namespace \ No newline at end of file diff --git a/interfaces/Matlab/cpp_classes/MxStruct.hpp b/interfaces/Matlab/cpp_classes/MxStruct.hpp new file mode 100644 index 000000000..924a81120 --- /dev/null +++ b/interfaces/Matlab/cpp_classes/MxStruct.hpp @@ -0,0 +1,29 @@ +// Copyright (c) 2025 Stefano Lovato and Charlie Vanaret +// Licensed under the MIT license. See LICENSE file in the project directory for details. + +#ifndef UNO_MXSTRUCT_H +#define UNO_MXSTRUCT_H + +#include +#include +#include "mex.h" + +namespace uno { + + class MxStruct { + + public: + mxArray* operator[](const std::string& key) const; + void insert(const std::string& key, mxArray* value); + [[nodiscard]] bool contains(const std::string& key) const; + [[nodiscard]] auto begin() const { return this->data.begin(); } + [[nodiscard]] auto end() const { return this->data.end(); } + + private: + std::map data; + }; + +}; // namespace + + +#endif // UNO_MXSTRUCT_H \ No newline at end of file diff --git a/interfaces/Matlab/example/example_hs015.m b/interfaces/Matlab/example/example_hs015.m new file mode 100644 index 000000000..32e02f253 --- /dev/null +++ b/interfaces/Matlab/example/example_hs015.m @@ -0,0 +1,132 @@ +%% HS015 problem (low-level interface) +clc, clear + +%% Optimization model +% Problem type: 'LP' = Linear, 'QP' = Quadratic, 'NLP' = Nonlinear +model.problem_type = 'NLP'; +% Vector indexing +model.base_indexing = 1; +% Number of variables +model.number_variables = 2; +% Variable bounds +model.variables_lower_bounds = [-inf, -inf]; +model.variables_upper_bounds = [0.5, inf]; + +% Optimization sense: 1 = Minimize, -1 = Maximize +model.optimization_sense = 1; +% Objective function handle: objective = objective_function(x) +model.objective_function = @objective_function; +% Objective gradient handle: gradient = objective_gradient(x) +model.objective_gradient = @objective_gradient; + +% Number of constraints +model.number_constraints = 2; +% Constraint bounds +model.constraints_lower_bounds = [1.0, 0.0]; +model.constraints_upper_bounds = [inf, inf]; +% Constraint function handle: constraints = constrain_function(x) +model.constraint_function = @constraint_function; +% Constraint jacobian handle: jacobian = constraint_jacobian(x) +model.constraint_jacobian = @constraint_jacobian; +% Constraint sparsity pattern (base_indexing-based) +model.number_jacobian_nonzeros = 4; +model.jacobian_row_indices = [1 2 1 2]; +model.jacobian_column_indices = [1 1 2 2]; + +% Lagrangian sign convention: 1 = rho*f(x) + y^T c(x), -1 = rho*f(x) - y^T c(x) +model.lagrangian_sign_convention = -1; +% Hessian triangular part: 'L' = lower, 'U' = upper +model.hessian_triangular_part = 'L'; +% Lagrangian Hessian handle: hessian = lagrangian_hessian(x,rho,y) +model.lagrangian_hessian = @lagrangian_hessian; +% Lagrangian Hessian sparsity pattern (base_indexing-based) +model.number_hessian_nonzeros = 3; +model.hessian_row_indices = [1 2 2]; +model.hessian_column_indices = [1 1 2]; + +% Jacobian operator handle: result = jacobian_operator(x,v) +% model.jacobian_operator = @jacobian_operator; + +% Jacobian tranposed operator handle: result = jacobian_transposed_operator(x,v) +% model.jacobian_transposed_operator = @jacobian_transposed_operator; + +% Hessian operator handle: result = lagrangian_hessian_operator(x,rho,y,v) +% model.lagrangian_hessian_operator = @lagrangian_hessian_operator; + +% Initial primal point +model.initial_primal_iterate = [-2.0, 1.0]; +% Dual primal point +model.initial_dual_iterate = [0, 0]; + +%% UNO solver options +% Preset can be possibly given: options = uno_options([preset]) +options = uno_options(); % default preset +% options = uno_options('ipopt'); + +%% UNO solver callbacks +callbacks = struct(); + +% Logger stream handle: logger_stream_callback(str) +% callbacks.logger_stream = @logger_stream_callback; + +% Notify acceptable iterate handle: notify_acceptable_iterate_callback(x, yl, yb, y, rho, feas, stat, compl) +% callbacks.notify_acceptable_iterate = @notify_acceptable_iterate_callback; + +% User termination handle: terminate = termination_callback(x, yl, yb, y, rho, feas, stat, compl) +% callbacks.termination = @termination_callback; + +%% Solving the model +result = uno_optimize(model, options, callbacks); +% Clear mex functions +clear uno_optimize uno_options + +%% Display the result +% disp(result) + +%% Model functions +% Objective function +function fval = objective_function(x) + fval = 100 * (x(2) - x(1)^2)^2 + (1 - x(1))^2; +end + +% Constraint function +function c = constraint_function(x) + c = zeros(2,1); + c(1) = x(1)*x(2); + c(2) = x(1) + x(2)^2; +end + +% Objective gradient +function grad = objective_gradient(x) + grad = zeros(2,1); + grad(1) = 400*x(1)^3 - 400*x(1)*x(2) + 2*x(1) - 2; + grad(2) = 200*(x(2) - x(1)^2); +end + +% Constraint jacobian +function J = constraint_jacobian(x) + J = zeros(4,1); + J(1) = x(2); + J(2) = 1; + J(3) = x(1); + J(4) = 2*x(2); +end + +% Lagrangian Hessian (lower triangular) +function H = lagrangian_hessian(x, rho, y) + H = zeros(3,1); + H(1) = rho*(1200*x(1)^2 - 400*x(2) + 2); + H(2) = -400*rho*x(1) - y(1); + H(3) = 200*rho - 2*y(2); +end + +% Lagrangian Hessian operator +function result = lagrangian_hessian_operator(x, rho, y, v) + h00 = rho*(1200*x(1)^2 - 400*x(2) + 2); + h10 = -400*rho*x(1) - y(1); + h11 = 200*rho - 2*y(2); + + result = zeros(2,1); + result(1) = h00*v(1) + h10*v(2); + result(2) = h10*v(1) + h11*v(2); +end \ No newline at end of file diff --git a/interfaces/Matlab/example/example_options.m b/interfaces/Matlab/example/example_options.m new file mode 100644 index 000000000..6c13de9fe --- /dev/null +++ b/interfaces/Matlab/example/example_options.m @@ -0,0 +1,8 @@ +%% UNO's option example +clc, clear + +% default preset +options_default = uno_options(); + +% IPOPT preset +options_ipopt = uno_options('ipopt'); \ No newline at end of file diff --git a/interfaces/Matlab/example/example_polak5.m b/interfaces/Matlab/example/example_polak5.m new file mode 100644 index 000000000..3b70790ab --- /dev/null +++ b/interfaces/Matlab/example/example_polak5.m @@ -0,0 +1,132 @@ +%% Polak5 prolbem (low-level interface) +clc, clear + +%% Optimization model +% Problem type: 'LP' = Linear, 'QP' = Quadratic, 'NLP' = Nonlinear +model.problem_type = 'NLP'; +% Vector indexing +model.base_indexing = 1; +% Number of variables +model.number_variables = 3; +% Variable bounds +model.variables_lower_bounds = [-inf, -inf, -inf]; +model.variables_upper_bounds = [ inf, inf, inf]; + +% Optimization sense: 1 = Minimize, -1 = Maximize +model.optimization_sense = 1; +% Objective function handle: objective = objective_function(x) +model.objective_function = @objective_function; +% Objective gradient handle: gradient = objective_gradient(x) +model.objective_gradient = @objective_gradient; + +% Number of constraints +model.number_constraints = 2; +% Constraint bounds +model.constraints_lower_bounds = [-inf, -inf]; +model.constraints_upper_bounds = [0, 0]; +% Constraint function handle: constraints = constrain_function(x) +model.constraint_function = @constraint_function; +% Constraint jacobian handle: jacobian = constraint_jacobian(x) +model.constraint_jacobian = @constraint_jacobian; +% Constraint sparsity pattern (base_indexing-based) +model.number_jacobian_nonzeros = 6; +model.jacobian_row_indices = [1 1 1 2 2 2]; +model.jacobian_column_indices = [1 2 3 1 2 3]; + +% Lagrangian sign convention: 1 = rho*f(x) + y^T c(x), -1 = rho*f(x) - y^T c(x) +model.lagrangian_sign_convention = -1; +% Hessian triangular part: 'L' = lower, 'U' = upper +model.hessian_triangular_part = 'L'; +% Lagrangian Hessian handle: hessian = lagrangian_hessian(x,rho,y) +model.lagrangian_hessian = @lagrangian_hessian; +% Lagrangian Hessian sparsity pattern (base_indexing-based) +model.number_hessian_nonzeros = 3; +model.hessian_row_indices = [1 2 2]; +model.hessian_column_indices = [1 1 2]; + +% Jacobian operator handle: result = jacobian_operator(x,v) +% model.jacobian_operator = @jacobian_operator; + +% Jacobian tranposed operator handle: result = jacobian_transposed_operator(x,v) +% model.jacobian_transposed_operator = @jacobian_transposed_operator; + +% Hessian operator handle: result = lagrangian_hessian_operator(x,rho,y,v) +% model.lagrangian_hessian_operator = @lagrangian_hessian_operator; + +% Initial primal point +model.initial_primal_iterate = [0.1, 0.1, 0.0]; +% Dual primal point +model.initial_dual_iterate = [0, 0]; + +%% UNO solver options +% Preset can be possibly given: options = uno_options([preset]) +options = uno_options(); % default preset +% options = uno_options('ipopt'); + +%% UNO solver callbacks +callbacks = struct(); + +% Logger stream handle: logger_stream_callback(str) +% callbacks.logger_stream = @logger_stream_callback; + +% Notify acceptable iterate handle: notify_acceptable_iterate_callback(x, yl, yb, y, rho, feas, stat, compl) +% callbacks.notify_acceptable_iterate = @notify_acceptable_iterate_callback; + +% User termination handle: terminate = termination_callback(x, yl, yb, y, rho, feas, stat, compl) +% callbacks.termination = @termination_callback; + +%% Solving the model +result = uno_optimize(model, options, callbacks); +% Clear mex functions +clear uno_optimize uno_options + +%% Display the result +% disp(result) + +%% Model functions +% Objective +function fval = objective_function(x) + fval = x(3); +end + +% Objective gradient +function grad = objective_gradient(~) + grad = [0; 0; 1]; +end + +% Constraint functions +function c = constraint_function(x) + c = zeros(2,1); + c(1) = -x(3) + 3*x(1)^2 + 50*(x(1) - x(2)^4 - 1)^2; + c(2) = -x(3) + 3*x(1)^2 + 50*(x(1) - x(2)^4 + 1)^2; +end + +% Constraint Jacobian +function J = constraint_jacobian(x) + J = zeros(6,1); + J(1) = 6*x(1) + 100*(x(1) - x(2)^4 - 1); + J(2) = -400*x(2)^3*(x(1) - x(2)^4 - 1); + J(3) = -1; + J(4) = 6*x(1) + 100*(x(1) - x(2)^4 + 1); + J(5) = -400*x(2)^3*(x(1) - x(2)^4 + 1); + J(6) = -1; +end + +% Lagrangian Hessian (lower triangular) +function H = lagrangian_hessian(x, ~, y) + H = zeros(3,1); + H(1) = -106*y(1) - 106*y(2); + H(2) = 400*x(2)^3*y(1) + 400*x(2)^3*y(2); + H(3) = 1200*x(2)^2*y(2)*(- x(2)^4 + x(1) + 1) - 1600*x(2)^6*y(2) - 1200*x(2)^2*y(1)*(x(2)^4 - x(1) + 1) - 1600*x(2)^6*y(1); +end + +% Lagrangian Hessian operator +function result = lagrangian_hessian_operator(x, ~, y, v) + h00 = -106*y(1) - 106*y(2); + h10 = 400*x(2)^3*y(1) + 400*x(2)^3*y(2); + h11 = 1200*x(2)^2*y(2)*(- x(2)^4 + x(1) + 1) - 1600*x(2)^6*y(2) - 1200*x(2)^2*y(1)*(x(2)^4 - x(1) + 1) - 1600*x(2)^6*y(1); + + result = zeros(3,1); + result(1) = h00*v(1) + h10*v(2); + result(2) = h10*v(1) + h11*v(2); +end \ No newline at end of file diff --git a/interfaces/Matlab/example/example_uno.m b/interfaces/Matlab/example/example_uno.m new file mode 100644 index 000000000..b3bed4b8b --- /dev/null +++ b/interfaces/Matlab/example/example_uno.m @@ -0,0 +1,96 @@ +%% UNO high-level interface example +clc, clear + +%% HS015 problem +lb = [-inf, -inf]; +ub = [0.5, inf]; +x0 = [-2.0, 1.0]; +% options +options = struct(); +options.HessianFnc = @hessian_hs015; +% call uno +[x,fval,exitflag,output,lambda,grad,hessian] = uno(@fun_hs015,x0,[],[],[],[],lb,ub,@nlcon_hs015,options); + +%% Polak5 problem +x0 = [0.1, 0.1, 0.0]; +% options +options = struct(); +% options.HessianFnc = @hessian_polak5; % hessian_model='LBFGS' used if HessianFnc is not provided +% call uno +[x,fval,exitflag,output,lambda,grad,hessian] = uno(@fun_polak5,x0,[],[],[],[],[],[],@nlcon_polak5,options); + +%% HS015 functions +% Objective function +% objective gradient must be also specified +function [fval,grad] = fun_hs015(x) + fval = 100 * (x(2) - x(1)^2)^2 + (1 - x(1))^2; + if nargout > 1 + grad(1) = 400*x(1)^3 - 400*x(1)*x(2) + 2*x(1) - 2; + grad(2) = 200*(x(2) - x(1)^2); + end +end + +% Constraint function +% constraint gradient must be also specified +function [c,ceq,gradc,gradceq] = nlcon_hs015(x) + c(1) = -x(1)*x(2)+1; + c(2) = -(x(1) + x(2)^2); + ceq = []; + if nargout > 2 + gradc = [-x(2), -1; ... + -x(1), -2*x(2)]; + gradceq = []; + end +end + +% Hessian function (optional) +function H = hessian_hs015(x,rho,lambda) + % objective + H = [1200*x(1)^2 - 400*x(2) + 2, -400*x(1); + -400*x(1), 200] * rho; + % c(1) + H = H + [0, -1; + -1, 0] * lambda.ineqnonlin(1); + % c(2) + H = H + [0, 0; + 0, 1] * lambda.ineqnonlin(2); +end + +%% Polak5 functions +% Objective function +% objective gradient and hessian must be also specified +function [fval,grad] = fun_polak5(x) + fval = x(3); + if nargout > 1 + grad = [0; 0; 1]; + end +end + +% Constraint function +% constraint gradient and hessian must be also specified +% constraint hessian is specified as a multidimensional matrix +function [c,ceq,gradc,gradceq] = nlcon_polak5(x) + c(1) = -x(3) + 3*x(1)^2 + 50*(x(1) - x(2)^4 - 1)^2; + c(2) = -x(3) + 3*x(1)^2 + 50*(x(1) - x(2)^4 + 1)^2; + ceq = []; + if nargout > 2 + gradc = [6*x(1) + 100*(x(1) - x(2)^4 - 1), 6*x(1) + 100*(x(1) - x(2)^4 + 1); ... + -400*x(2)^3*(x(1) - x(2)^4 - 1), -400*x(2)^3*(x(1) - x(2)^4 + 1); ... + -1, -1]; + gradceq = []; + end +end + +% Hessian function (optional) +function H = hessian_polak5(x,rho,lambda) + % objective + H = zeros(3,3) * rho; + % c(1) + H = H + [106, -400*x(2)^3, 0; + -400*x(2)^3, -1200*x(2)^2*(x(1)-x(2)^4-1)+1600*x(2)^6, 0; + 0, 0, 0] * lambda.ineqnonlin(1); + % c(2) + H = H + [106, -400*x(2)^3, 0; ... + -400*x(2)^3, -1200*x(2)^2*(x(1)-x(2)^4+1)+1600*x(2)^6, 0; + 0, 0, 0] * lambda.ineqnonlin(2); +end \ No newline at end of file diff --git a/interfaces/Matlab/uno.m b/interfaces/Matlab/uno.m new file mode 100644 index 000000000..2dca059f6 --- /dev/null +++ b/interfaces/Matlab/uno.m @@ -0,0 +1,360 @@ +function [X,FVAL,EXITFLAG,OUTPUT,LAMBDA,GRAD,HESSIAN] = uno(FUN,X,A,B,Aeq,Beq,LB,UB,NONLCON,options,varargin) +% UNO - UNO high-level interface +% Solve an optimization model described using UNO, a modern and modular solver for nonlinearly constrained optimization. +% +% Syntax +% x = uno(fun,x0) +% x = uno(fun,x0,A,b) +% x = uno(fun,x0,A,b,Aeq,beq) +% x = uno(fun,x0,A,b,Aeq,beq,lb,ub) +% x = uno(fun,x0,A,b,Aeq,beq,lb,ub,nonlcon) +% x = uno(fun,x0,A,b,Aeq,beq,lb,ub,nonlcon,options) +% options = uno('defaults') +% options = uno('defaults', preset) +% +% Input Arguments +% fun - Function to minimize +% function handle | function name +% x0 - Initial point +% real vector | real array +% A - Linear inequality constraints +% real matrix +% b - Linear inequality constraints +% real vector +% Aeq - Linear equality constraints +% real matrix +% beq - Linear equality constraints +% real vector +% lb - Lower bounds +% real vector | real array +% ub - Upper bounds +% real vector | real array +% nonlcon - Nonlinear constraints +% function handle | function name +% options - UNO options, as returned by uno('defaults') +% structure +% preset - Option preset +% character vector | string +% +% Output Arguments +% x - Solution +% real vector | real array +% fval - Objective function value at solution +% real number +% exitflag - Reason fmincon stopped +% integer +% output - Information about the optimization process +% structure +% lambda - Lagrange multipliers at the solution +% structure +% grad - Gradient at the solution +% real vector +% hessian - Hessian matrix +% real matrix +% options - UNO options +% structure +% +% Examples +% example_uno +% +% Copyright (c) 2025 Stefano Lovato and Charlie Vanaret +% Licensed under the MIT license. See LICENSE file in the project directory for details. + +% return the options if 'defaults' and preset +if nargin<=2 && nargout <= 1 && strcmpi(FUN,'defaults') + if nargin==1 + X = uno_options(); + else + X = uno_options(X); + end + return +end + +% default args +if nargin < 10 + options = struct(); + if nargin < 9 + NONLCON = []; + if nargin < 8 + UB = []; + if nargin < 7 + LB = []; + if nargin < 6 + Beq = []; + if nargin < 5 + Aeq = []; + if nargin < 4 + B = []; + if nargin < 3 + A = []; + end + end + end + end + end + end + end +end + +% Column vector +X = X(:); +Beq = Beq(:); +B = B(:); + +% Check for complex X0 +if ~isreal(X) + error('Initial point contains complex values. Uno cannot continue.') +end + +% Check for empty X +if isempty(X) + error('Initial point must be non-empty.'); +end + +% Set empty linear constraints +if isempty(Aeq) + Aeq = zeros(0,length(X)); +end +if isempty(A) + A = zeros(0,length(X)); +end +if isempty(Beq) + Beq = zeros(0,1); +end +if isempty(B) + B = zeros(0,1); +end + +% Check linear constraint sizes +if size(Aeq,2) ~= length(X) + error('Aeq must have %d column(s)', length(X)); +end +if size(Aeq,1) ~= length(Beq) + error('Row dimension of Aeq is inconsistent with length of beq.'); +end +if size(A,2) ~= length(X) + error('A must have %d column(s)', length(X)); +end +if size(A,1) ~= length(B) + error('Row dimension of A is inconsistent with length of b.'); +end + +% Check the bounds +LB = LB(:); +UB = UB(:); +if length(LB) > length(X) + LB = LB(1:length(X)); + warning('Length of lower bounds is > length(x); ignoring extra bounds.') +end +if length(UB) > length(X) + UB = UB(1:length(X)); + warning('Length of upper bounds is > length(x); ignoring extra bounds.') +end +if any (UB == -Inf) + error('-Inf detected in upper bound: upper bounds must be > -Inf.'); +end +if any (LB == +Inf) + error('+Inf detected in upper bound: upper bounds must be < Inf.'); +end +if anynan(LB) || anynan(UB) + error('Lower and upper bounds must not be NaN.'); +end +LB = [LB; -inf(length(X)-length(LB),1)]; +UB = [UB; +inf(length(X)-length(UB),1)]; + +% Check [fval,grad] = FUN(x) +if ~isa(FUN,'function_handle') + error('FUN must be a function handle.'); +end +if (nargin(FUN)~=1) || (nargout(FUN)<2) + error('FUN must have 1 input and 2 output argument(s).') +end +try + [fval0,grad0] = FUN(X); +catch + error('Failure in objective function evaluation.') +end +if ~isnumeric(fval0) || ~isnumeric(grad0) + error('FUM must return numeric values.') +end +if ~isscalar(fval0) + error('FUM must return a scalar value.') +end +if length(grad0) ~= length(X) + error('Dimension of grad is inconsistent with the length of x.') +end + +% Check [c,ceq,gradc,gradceq] = NONLCON(x) +if isempty(NONLCON) + NONLCON = @(x) deal([]); +else + if ~isa(NONLCON,'function_handle') + error('NLCON must be a function handle.') + end + if (nargin(NONLCON)~=1) || (nargout(NONLCON)<4) + error('NONLCON must have 1 input and 4 output argument(s).') + end +end +try + [c0,ceq0,gradc0,gradceq0] = NONLCON(X); +catch + error('Failure in nonlinear constraint function evaluation.') +end +if ~isnumeric(c0) || ~isnumeric(ceq0) || ~isnumeric(gradc0) || ~isnumeric(gradceq0) + error('NONLCON must return numeric values.') +end +if length(c0) ~= size(gradc0,2) + error('Column dimension of gradc is inconsistent with the length of c.') +end +if length(ceq0) ~= size(gradceq0,2) + error('Column dimension of gradceq is inconsistent with the length of ceq.') +end +if ~isempty(c0) && length(X) ~= size(gradc0,1) + error('Row dimension of gradc is inconsistent with the length of x.') +end +if ~isempty(ceq0) && length(X) ~= size(gradceq0,1) + error('Row dimension of gradceq is inconsistent with the length of x.') +end + +% Check H = HESSIANFNC(x,rho,lambda) +HESSIANFNC = []; +if isfield(options,'HessianFnc') + HESSIANFNC = options.HessianFnc; + options = rmfield(options, 'HessianFnc'); + if ~isa(HESSIANFNC,'function_handle') + error('HESSIANFNC must be a function handle.') + end + if (nargin(HESSIANFNC)~=3) || (nargout(HESSIANFNC)<1) + error('HESSIANFNC must have 3 input and 1 output argument(s).') + end + LAMBDA.eqnonlin = zeros(length(ceq0),1); + LAMBDA.ineqnonlin = zeros(length(c0),1); + try + H0 = HESSIANFNC(X,1,LAMBDA); + catch + error('Failure in Lagrangian Hessian function evaluation.') + end + if ~isnumeric(H0) + error('HESSIANFNC must return numeric values.') + end + if ~isempty(H0) && ~isequal(size(H0), [length(X), length(X)]) + error('Dimension of Lagrangian Hessian is inconsistent with the length of x.') + end +end + +% Constraint indexes +ind_eqlin = 1 : numel(Beq); +ind_ineqlin = (numel(Beq)+1) : (numel(Beq)+numel(B)); +ind_ineqnonlin = (numel(Beq)+numel(B)+1) : (numel(Beq)+numel(B)+numel(c0)); +ind_eqnonlin = (numel(Beq)+numel(B)+numel(c0)+1) : (numel(Beq)+numel(B)+numel(c0)+numel(ceq0)); + +% Create the optimization problem +model.problem_type = 'NLP'; +model.base_indexing = 1; +model.number_variables = length(X); +model.variables_lower_bounds = LB; +model.variables_upper_bounds = UB; +model.optimization_sense = 1; +model.objective_function = @(x) objective_function(x, FUN); +model.objective_gradient = @(x) objective_gradient(x, FUN); +model.number_constraints = length(Beq) + length(B) + length(c0) + length(ceq0); +model.constraints_lower_bounds = [zeros(size(Beq)); -inf(size(B)); -inf(length(c0),1); zeros(length(ceq0),1)]; +model.constraints_upper_bounds = zeros(model.number_constraints,1); +model.constraint_function = @(x) constraint_function(x,Aeq,Beq,A,B,NONLCON); +[irj, jcj] = find(ones(model.number_constraints,length(X))); +model.jacobian_row_indices = irj; +model.jacobian_column_indices = jcj; +model.number_jacobian_nonzeros = length(irj); +model.constraint_jacobian = @(x) constraint_jacobian(x,Aeq,Beq,A,B,NONLCON,irj,jcj); +model.lagrangian_sign_convention = 1; +model.initial_primal_iterate = X; + +% Lagrangian Hessian +if ~isempty(HESSIANFNC) + model.hessian_triangular_part = 'L'; + [irh, jch] = find(tril(ones(length(X),length(X)))); + model.hessian_row_indices = irh(:); + model.hessian_column_indices = jch(:); + model.number_hessian_nonzeros = length(irh(:)); + model.lagrangian_hessian = @(x,rho,y) lagrangian_hessian(x,rho,y,HESSIANFNC,ind_ineqnonlin,ind_eqnonlin,irh,jch); +else + if ~isfield(options,'hessian_model') || strcmp(options.hessian_model,'exact') + % fallback to LBFGS when no HessianFnc + options.hessian_model = 'LBFGS'; + end +end + +% Call UNO +result = uno_optimize(model, options); + +% Results +X = result.primal_solution(1:length(X)); +FVAL = result.solution_objective; +switch result.optimization_status + case 0 + switch result.solution_status + case {0,3,5,6} + EXITFLAG = -2; + case {1,2} + EXITFLAG = 1; + case 4 + EXITFLAG = 2; + end + case 1 + EXITFLAG = 0; + case {2,3,4} + EXITFLAG = -2; + case 5 + EXITFLAG = -1; +end +OUTPUT.iterations = result.number_iterations; +OUTPUT.funcCount = result.number_objective_evaluations; +OUTPUT.constrviolation = result.solution_primal_feasibility; +OUTPUT.firstorderopt = result.solution_stationarity; +LAMBDA.lower = result.lower_bound_dual_solution(1:length(X)); +LAMBDA.upper = result.upper_bound_dual_solution(1:length(X)); +LAMBDA.eqlin = result.constraint_dual_solution(ind_eqlin); +LAMBDA.ineqlin = result.constraint_dual_solution(ind_ineqlin); +LAMBDA.ineqnonlin = result.constraint_dual_solution(ind_ineqnonlin); +LAMBDA.eqnonlin = result.constraint_dual_solution(ind_eqnonlin); +[~,GRAD] = FUN(X); +if ~isempty(HESSIANFNC) + HESSIAN = HESSIANFNC(X, 1, LAMBDA); +else + HESSIAN = []; +end + +end + +% Model functions +function fval = objective_function(x,FUN) + fval = FUN(x); +end + +function grad = objective_gradient(x,FUN) + [~, grad] = FUN(x); +end + +function constr = constraint_function(x,Aeq,Beq,A,B,NONLCON) + [c,ceq] = NONLCON(x); + constr = [Aeq*x-Beq; ... + A*x-B; ... + c(:); ... + ceq(:)]; +end + +function J = constraint_jacobian(x,Aeq,~,A,~,NONLCON,ir,jc) + [~,~,gradc,gradceq] = NONLCON(x); + J = [Aeq; ... + A; ... + gradc'; ... + gradceq']; + J = J(ir + (jc-1)*size(J,1)); +end + +function H = lagrangian_hessian(x, rho, y, HESSIANFNC, icneq, iceq, ir, jc) + LAMBDA.ineqlin = y(iceq); + LAMBDA.ineqnonlin = y(icneq); + H = HESSIANFNC(x, rho, LAMBDA); + H = H(ir + (jc-1)*size(H,1)); +end \ No newline at end of file diff --git a/interfaces/Matlab/uno_optimize.cpp b/interfaces/Matlab/uno_optimize.cpp new file mode 100644 index 000000000..745f82ed0 --- /dev/null +++ b/interfaces/Matlab/uno_optimize.cpp @@ -0,0 +1,343 @@ +// Copyright (c) 2025 Stefano Lovato and Charlie Vanaret +// Licensed under the MIT license. See LICENSE file in the project directory for details. + +#include +#include "Uno.hpp" +#include "options/DefaultOptions.hpp" +#include "options/Presets.hpp" +#include "tools/Logger.hpp" +#include "cpp_classes/MatlabModel.hpp" +#include "cpp_classes/MatlabUserCallbacks.hpp" +#include "cpp_classes/MatlabStreamCallback.hpp" +#include "cpp_classes/MxStruct.hpp" +#include "cpp_classes/ErrorString.hpp" +#include "unomex/unomex_conversion.hpp" +#include "unomex/unomex_validation.hpp" +#include "unomex/unomex_utils.hpp" +#include "mex.h" + +using namespace uno; + +// gateway function: result = uno_optimize(model[, options, callbacks]) +void mexFunction( int /* nlhs */, mxArray* plhs[], int nrhs, const mxArray* prhs[]) { +#if defined(__APPLE__) + // lock the mex to prevent matlab crash during cleanup + // cf. https://github.com/wmutschl/matlab_openmp_tls_bug + static bool locked = false; + if (!locked) + { + mexLock(); + locked = true; + } +#endif + + std::string errmsg; // declare once for all + // validate argument list + if (nrhs < 1) { + errmsg = ErrorString::format_error(ErrorType::NARGIN_NOTENOUGH, 1); + mexErrMsgIdAndTxt("uno:error", errmsg.c_str()); + } + if (nrhs > 3) { + errmsg = ErrorString::format_error(ErrorType::NARGIN_TOOMANY); + mexErrMsgIdAndTxt("uno:error", errmsg.c_str()); + } + + // model (mandatory) + if (!validate_struct_input(prhs[0], 1, errmsg)) { + mexErrMsgIdAndTxt("uno:error", errmsg.c_str()); + } + const MxStruct model = mxArray_to_mxStruct(prhs[0]); + + // options (optional) + MxStruct options; + if (nrhs > 1) { + if (!validate_struct_input(prhs[1], 2, errmsg)) { + mexErrMsgIdAndTxt("uno:error", errmsg.c_str()); + } + options = mxArray_to_mxStruct(prhs[1]); + } + + // callbacks (optional) + MxStruct callbacks; + if (nrhs > 2) { + if (!validate_struct_input(prhs[2], 3, errmsg)) { + mexErrMsgIdAndTxt("uno:error", errmsg.c_str()); + } + callbacks = mxArray_to_mxStruct(prhs[2]); + } + + // problem type (mandatory) + mxArray* problem_type = model["problem_type"]; + if (!validate_string_field(problem_type, {UNO_PROBLEM_LINEAR, UNO_PROBLEM_QUADRATIC, UNO_PROBLEM_NONLINEAR}, "problem_type", errmsg)) { + mexErrMsgIdAndTxt("uno:error", "Invalid model. %s", errmsg.c_str()); + } + + // base indexing (mandatory) + mxArray* base_indexing = model["base_indexing"]; + if (!validate_positive_integer_field(base_indexing, "base_indexing", errmsg)) { + mexErrMsgIdAndTxt("uno:error", "Invalid model. %s", errmsg.c_str()); + } + + // variables (mandatory) + mxArray* number_variables = model["number_variables"]; + mxArray* variables_lower_bounds = model["variables_lower_bounds"]; + mxArray* variables_upper_bounds = model["variables_upper_bounds"]; + if (!validate_positive_integer_field(number_variables, "number_variables", errmsg)) { + mexErrMsgIdAndTxt("uno:error", "Invalid model. %s", errmsg.c_str()); + } + if (!validate_double_vector_field(variables_lower_bounds, static_cast(mxArray_to_scalar(number_variables)), "variables_lower_bounds", errmsg)) { + mexErrMsgIdAndTxt("uno:error", "Invalid model. %s", errmsg.c_str()); + } + if (!validate_double_vector_field(variables_upper_bounds, static_cast(mxArray_to_scalar(number_variables)), "variables_lower_bounds", errmsg)) { + mexErrMsgIdAndTxt("uno:error", "Invalid model. %s", errmsg.c_str()); + } + + // objective (mandatory) + mxArray* objective_function = model["objective_function"]; + mxArray* objective_gradient = model["objective_gradient"]; + mxArray* optimization_sense = model["optimization_sense"]; + // objective = objective_function(x) + if (!validate_matlab_handle_field(objective_function, "objective_function", errmsg)) { + mexErrMsgIdAndTxt("uno:error", "Invalid model objective function. %s", errmsg.c_str()); + } + // gradient = objective_gradient(x) + if (!validate_matlab_handle_field(objective_gradient, "objective_gradient", errmsg)) { + mexErrMsgIdAndTxt("uno:error", "Invalid model objective gradient. %s", errmsg.c_str()); + } + if (!validate_unitary_field(optimization_sense, "optimization_sense", errmsg)) { + mexErrMsgIdAndTxt("uno:error", "Invalid model. %s", errmsg.c_str()); + } + + // constraints (mandatory) + mxArray* number_constraints = model["number_constraints"]; + mxArray* constraints_lower_bounds = model["constraints_lower_bounds"]; + mxArray* constraints_upper_bounds = model["constraints_upper_bounds"]; + mxArray* constraint_function = model["constraint_function"]; + if (!validate_positive_integer_field(number_constraints, "number_constraints", errmsg)) { + mexErrMsgIdAndTxt("uno:error", "Invalid model. %s", errmsg.c_str()); + } + if (!validate_double_vector_field(constraints_lower_bounds, static_cast(mxArray_to_scalar(number_constraints)), "constraints_lower_bounds", errmsg)) { + mexErrMsgIdAndTxt("uno:error", "Invalid model. %s", errmsg.c_str()); + } + if (!validate_double_vector_field(constraints_upper_bounds, static_cast(mxArray_to_scalar(number_constraints)), "constraints_upper_bounds", errmsg)) { + mexErrMsgIdAndTxt("uno:error", "Invalid model. %s", errmsg.c_str()); + } + // constraints = constrain_function(x) + if (!validate_matlab_handle_field(constraint_function, "constraint_function", errmsg)) { + mexErrMsgIdAndTxt("uno:error", "Invalid model constraint function. %s", errmsg.c_str()); + } + + // jacobian (mandatory) + mxArray* number_jacobian_nonzeros = model["number_jacobian_nonzeros"]; + mxArray* jacobian_row_indices = model["jacobian_row_indices"]; + mxArray* jacobian_column_indices = model["jacobian_column_indices"]; + mxArray* constraint_jacobian = model["constraint_jacobian"]; + if (!validate_positive_integer_field(number_jacobian_nonzeros, "number_jacobian_nonzeros", errmsg)) { + mexErrMsgIdAndTxt("uno:error", "Invalid model. %s", errmsg.c_str()); + } + if (!validate_integer_vector_field(jacobian_row_indices, static_cast(mxArray_to_scalar(number_jacobian_nonzeros)), "jacobian_row_indices", errmsg)) { + mexErrMsgIdAndTxt("uno:error", "Invalid model. %s", errmsg.c_str()); + } + if (!validate_integer_vector_field(jacobian_column_indices, static_cast(mxArray_to_scalar(number_jacobian_nonzeros)), "jacobian_column_indices", errmsg)) { + mexErrMsgIdAndTxt("uno:error", "Invalid model. %s", errmsg.c_str()); + } + // jacobian = constraint_jacobian(x) + if (!validate_matlab_handle_field(constraint_jacobian, "constraint_jacobian", errmsg)) { + mexErrMsgIdAndTxt("uno:error", "Invalid model constraint Jacobian. %s", errmsg.c_str()); + } + + // jacobian operators (optional) + mxArray* jacobian_operator = model["jacobian_operator"]; + mxArray* jacobian_transposed_operator = model["jacobian_transposed_operator"]; + if (jacobian_operator != nullptr) { + // result = jacobian_operator(x,v) + if (!validate_matlab_handle_field(jacobian_operator, "jacobian_operator", errmsg)) { + mexErrMsgIdAndTxt("uno:error", "Invalid model Jacobian operator. %s", errmsg.c_str()); + } + } + if (jacobian_transposed_operator != nullptr) { + // jacobian_transposed_operator(x,v) + if (!validate_matlab_handle_field(jacobian_transposed_operator, "jacobian_transposed_operator", errmsg)) { + mexErrMsgIdAndTxt("uno:error", "Invalid model Jacobian transposed operator. %s", errmsg.c_str()); + } + } + + // lagrangian sign convention (mandatory) + mxArray* lagrangian_sign_convention = model["lagrangian_sign_convention"]; + if (!validate_unitary_field(lagrangian_sign_convention, "lagrangian_sign_convention", errmsg)) { + mexErrMsgIdAndTxt("uno:error", "Invalid model. %s", errmsg.c_str()); + } + + // hessian (optional) + mxArray* hessian_triangular_part = model["hessian_triangular_part"]; + mxArray* number_hessian_nonzeros = model["number_hessian_nonzeros"]; + mxArray* hessian_row_indices = model["hessian_row_indices"]; + mxArray* hessian_column_indices = model["hessian_column_indices"]; + mxArray* lagrangian_hessian = model["lagrangian_hessian"]; + if (lagrangian_hessian != nullptr) { + if (!validate_char_field(hessian_triangular_part, {UNO_LOWER_TRIANGLE, UNO_UPPER_TRIANGLE}, "hessian_triangular_part", errmsg)) { + mexErrMsgIdAndTxt("uno:error", "Invalid model. %s", errmsg.c_str()); + } + if (!validate_positive_integer_field(number_hessian_nonzeros, "number_hessian_nonzeros", errmsg)) { + mexErrMsgIdAndTxt("uno:error", "Invalid model. %s", errmsg.c_str()); + } + if (!validate_integer_vector_field(hessian_row_indices, static_cast(mxArray_to_scalar(number_hessian_nonzeros)), "hessian_row_indices", errmsg)) { + mexErrMsgIdAndTxt("uno:error", "Invalid model. %s", errmsg.c_str()); + } + if (!validate_integer_vector_field(hessian_column_indices, static_cast(mxArray_to_scalar(number_hessian_nonzeros)), "hessian_column_indices", errmsg)) { + mexErrMsgIdAndTxt("uno:error", "Invalid model. %s", errmsg.c_str()); + } + // hessian = lagrangian_hessian(x,rho,y) + if (!validate_matlab_handle_field(lagrangian_hessian, "lagrangian_hessian", errmsg)) { + mexErrMsgIdAndTxt("uno:error", "Invalid model Lagrangian Hessian. %s", errmsg.c_str()); + } + } + + // hessian operators (optional) + mxArray* lagrangian_hessian_operator = model["lagrangian_hessian_operator"]; + if (lagrangian_hessian_operator != nullptr) { + // lagrangian_hessian_operator(x,rho,y,v) + if (!validate_matlab_handle_field(lagrangian_hessian_operator, "lagrangian_hessian_operator", errmsg)) { + mexErrMsgIdAndTxt("uno:error", "Invalid model Lagrangian Hessian operator. %s", errmsg.c_str()); + } + } + + // initial iterates (optional) + mxArray* initial_primal_iterate = model["initial_primal_iterate"]; + mxArray* initial_dual_iterate = model["initial_dual_iterate"]; + if (initial_primal_iterate != nullptr) { + if (!validate_double_vector_field(initial_primal_iterate, static_cast(mxArray_to_scalar(number_variables)), "initial_primal_iterate", errmsg)) { + mexErrMsgIdAndTxt("uno:error", "Invalid model. %s", errmsg.c_str()); + } + } + if (initial_dual_iterate != nullptr) { + if (!validate_double_vector_field(initial_dual_iterate, static_cast(mxArray_to_scalar(number_constraints)), "initial_dual_iterate", errmsg)) { + mexErrMsgIdAndTxt("uno:error", "Invalid model. %s", errmsg.c_str()); + } + } + + // callbacks (optional) + mxArray* logger_stream_callback = callbacks["logger_stream"]; + mxArray* notify_acceptable_iterate_callback = callbacks["notify_acceptable_iterate"]; + mxArray* termination_callback = callbacks["termination"]; + if (logger_stream_callback != nullptr) { + // logger_stream(str) + if (!validate_matlab_handle_field(logger_stream_callback, "logger_stream", errmsg)) { + mexErrMsgIdAndTxt("uno:error", "Invalid logger stream callback. %s", errmsg.c_str()); + } + } + if (notify_acceptable_iterate_callback != nullptr) { + // notify_acceptable_iterate(x, yl, yb, y, rho, feas, stat, compl) + if (!validate_matlab_handle_field(notify_acceptable_iterate_callback, "notify_acceptable_iterate", errmsg)) { + mexErrMsgIdAndTxt("uno:error", "Invalid notify acceptable iterate callback. %s", errmsg.c_str()); + } + } + if (termination_callback != nullptr) { + // terminate = termination(x, yl, yb, y, rho, feas, stat, compl) + if (!validate_matlab_handle_field(termination_callback, "termination", errmsg)) { + mexErrMsgIdAndTxt("uno:error", "Invalid user termination callback. %s", errmsg.c_str()); + } + } + + // set the logger stream + MatlabStreamCallback matlab_stream_callback(logger_stream_callback); + UserOStream matlab_ostream(&matlab_stream_callback); + Logger::set_stream(matlab_ostream); + + // create user callbacks + MatlabUserCallbacks user_callbacks(notify_acceptable_iterate_callback, termination_callback); + + // create user model + MatlabUserModel user_model(mxArray_to_string(problem_type).c_str(), + static_cast(mxArray_to_scalar(number_variables)), + static_cast(mxArray_to_scalar(base_indexing))); + user_model.number_constraints = static_cast(mxArray_to_scalar(number_constraints)); + user_model.variables_lower_bounds = mxArray_to_stdvector(variables_lower_bounds); + user_model.variables_upper_bounds = mxArray_to_stdvector(variables_upper_bounds); + // objective + user_model.objective_function = objective_function; + user_model.objective_gradient = objective_gradient; + user_model.optimization_sense = static_cast(mxArray_to_scalar(optimization_sense)); + // constraint + user_model.constraint_functions = constraint_function; + user_model.constraints_lower_bounds = mxArray_to_stdvector(constraints_lower_bounds); + user_model.constraints_upper_bounds = mxArray_to_stdvector(constraints_upper_bounds); + user_model.jacobian = constraint_jacobian; + user_model.number_jacobian_nonzeros = static_cast(mxArray_to_scalar(number_jacobian_nonzeros)); + Vector jacobian_row_indices_vector = convert_vector_type( mxArray_to_vector(jacobian_row_indices) ); + Vector jacobian_column_indices_vector = convert_vector_type( mxArray_to_vector(jacobian_column_indices) ); + user_model.jacobian_row_indices = std::vector(jacobian_row_indices_vector.begin(), jacobian_row_indices_vector.end()); + user_model.jacobian_column_indices = std::vector(jacobian_column_indices_vector.begin(), jacobian_column_indices_vector.end()); + // jacobian operators + user_model.jacobian_operator = jacobian_operator; + user_model.jacobian_transposed_operator = jacobian_transposed_operator; + // hessian + user_model.lagrangian_sign_convention = static_cast(mxArray_to_scalar(lagrangian_sign_convention)); + user_model.lagrangian_hessian = lagrangian_hessian; + if (user_model.lagrangian_hessian) { + user_model.hessian_triangular_part = mxArray_to_scalar(hessian_triangular_part); + user_model.number_hessian_nonzeros = static_cast(mxArray_to_scalar(number_hessian_nonzeros)); + Vector hessian_row_indices_vector = convert_vector_type( mxArray_to_vector(hessian_row_indices) ); + Vector hessian_column_indices_vector = convert_vector_type( mxArray_to_vector(hessian_column_indices) ); + user_model.hessian_row_indices = std::vector(hessian_row_indices_vector.begin(), hessian_row_indices_vector.end()); + user_model.hessian_column_indices = std::vector(hessian_column_indices_vector.begin(), hessian_column_indices_vector.end()); + } + // hessian operator + user_model.lagrangian_hessian_operator = lagrangian_hessian_operator; + // initial iterates + if (initial_primal_iterate != nullptr) { + user_model.initial_primal_iterate = mxArray_to_stdvector(initial_primal_iterate); + } else { + user_model.initial_primal_iterate = std::vector(user_model.number_variables, 0.); + } + if (initial_dual_iterate != nullptr) { + user_model.initial_dual_iterate = mxArray_to_stdvector(initial_dual_iterate); + } else { + user_model.initial_dual_iterate = std::vector(user_model.number_constraints, 0.); + } + + // create Uno model + const MatlabModel uno_model(user_model); + + // create Uno solver + Uno uno_solver; + + // create Uno options + Options uno_options; + DefaultOptions::load(uno_options); + // handle preset + Presets::set_default(uno_options); + if (options.contains("preset")) { + const std::string preset = mxArray_to_string(options["preset"]); + // set only if it is valid (not empty) + if (!preset.empty()) { + Presets::set(uno_options, preset); + } + } + // override with user options + mxStruct_to_options(options, uno_options); + + // check lagrangian hessian + if (uno_options.get_string("hessian_model") == std::string("exact") && + user_model.lagrangian_hessian == nullptr && user_model.lagrangian_hessian_operator == nullptr) { + mexErrMsgIdAndTxt("uno:error", "Lagrangian Hessian must be provided when hessian_model is 'exact'"); + } + + // solve + Logger::set_logger(uno_options.get_string("logger")); + MxStruct result; + try { + Result uno_result = uno_solver.solve(uno_model, uno_options, user_callbacks); + result = result_to_mxStruct(uno_result); + } + catch (const std::exception& e) { + errmsg = ErrorString::format_error(ErrorType::UNO, e.what()); + mexErrMsgIdAndTxt("uno:error", errmsg.c_str()); + } + + // output + plhs[0] = mxStruct_to_mxArray(result); + + // flush the logger + Logger::flush(); +} diff --git a/interfaces/Matlab/uno_optimize.m b/interfaces/Matlab/uno_optimize.m new file mode 100644 index 000000000..b5266c731 --- /dev/null +++ b/interfaces/Matlab/uno_optimize.m @@ -0,0 +1,29 @@ +% UNO_OPTIMIZE - UNO solver +% Call the UNO solver and solve the specified optimization model. +% +% Syntax +% result = UNO_OPTIMIZE(model) +% result = UNO_OPTIMIZE(model,options) +% result = UNO_OPTIMIZE(model,options,callbacks) +% +% Input Arguments +% model - Optimization model +% structure +% options - UNO options, as returned by uno_options +% structure +% callbacks - UNO solver callbacks +% structure +% +% Output Arguments +% result - UNO result +% structure +% +% Examples +% example_options +% example_hs015 +% example_polak5 +% +% See also uno_options uno +% +% Copyright (c) 2025 Stefano Lovato and Charlie Vanaret +% Licensed under the MIT license. See LICENSE file in the project directory for details. \ No newline at end of file diff --git a/interfaces/Matlab/uno_options.cpp b/interfaces/Matlab/uno_options.cpp new file mode 100644 index 000000000..d601cbea8 --- /dev/null +++ b/interfaces/Matlab/uno_options.cpp @@ -0,0 +1,62 @@ +// Copyright (c) 2025 Stefano Lovato and Charlie Vanaret +// Licensed under the MIT license. See LICENSE file in the project directory for details. + +#include "options/DefaultOptions.hpp" +#include "options/Presets.hpp" +#include "cpp_classes/ErrorString.hpp" +#include "unomex/unomex_conversion.hpp" +#include "unomex/unomex_utils.hpp" +#include "mex.h" + +using namespace uno; + +// gateway function: options = uno_options([preset]); +void mexFunction( int /* nlhs */, mxArray* plhs[], int nrhs, const mxArray* prhs[]) { +#if defined(__APPLE__) + // lock the mex to prevent matlab crash during cleanup + // cf. https://github.com/wmutschl/matlab_openmp_tls_bug + static bool locked = false; + if (!locked) + { + mexLock(); + locked = true; + } +#endif + + // validate arguments + if (nrhs > 1) { + const std::string errmsg = ErrorString::format_error(ErrorType::NARGIN_TOOMANY); + mexErrMsgIdAndTxt("uno:error", errmsg.c_str()); + } + // preset (optional) + std::string preset; + if (nrhs == 1 && !isempty(prhs[0])) { + if (!isa(prhs[0]) && !isa(prhs[0])) { + const std::string errmsg = ErrorString::format_error(ErrorType::INPUT_STRING, 1); + mexErrMsgIdAndTxt("uno:error", errmsg.c_str()); + } + preset = mxArray_to_string(prhs[0]); + } + + // create Uno options + Options uno_options; + DefaultOptions::load(uno_options); + // add default preset + Presets::set_default(uno_options); + // set preset + if (!preset.empty()) { + try { + Presets::set(uno_options, preset); + } + catch (const std::runtime_error& err) { + const std::string errmsg = ErrorString::format_error(ErrorType::INPUT_INVALID, 1); + mexErrMsgIdAndTxt("uno:error", "%s %s.", errmsg.c_str(), err.what()); + } + } + + // options + MxStruct options = options_to_mxStruct(uno_options); + + // output + plhs[0] = mxStruct_to_mxArray(options); +} \ No newline at end of file diff --git a/interfaces/Matlab/uno_options.m b/interfaces/Matlab/uno_options.m new file mode 100644 index 000000000..f8b94c68f --- /dev/null +++ b/interfaces/Matlab/uno_options.m @@ -0,0 +1,23 @@ +% UNO_OPTIONS - UNO options +% Get the UNO default options with optionally a given option preset. +% +% Syntax +% options = UNO_OPTIONS() +% options = UNO_OPTIONS(preset) +% +% Input Arguments +% preset - Option preset +% character vector | string +% +% Output Arguments +% options - UNO options +% structure +% +% Examples +% example_hs015 +% example_polak5 +% +% See also uno_optimize uno +% +% Copyright (c) 2025 Stefano Lovato and Charlie Vanaret +% Licensed under the MIT license. See LICENSE file in the project directory for details. \ No newline at end of file diff --git a/interfaces/Matlab/unomex/unomex_conversion.cpp b/interfaces/Matlab/unomex/unomex_conversion.cpp new file mode 100644 index 000000000..489d45a93 --- /dev/null +++ b/interfaces/Matlab/unomex/unomex_conversion.cpp @@ -0,0 +1,184 @@ +// Copyright (c) 2025 Stefano Lovato and Charlie Vanaret +// Licensed under the MIT license. See LICENSE file in the project directory for details. + +#include "tools/Logger.hpp" +#include "../cpp_classes/ErrorString.hpp" +#include "unomex_utils.hpp" +#include "unomex_conversion.hpp" + +namespace uno { + + std::string mxArray_to_string(const mxArray* arr) { + std::string str; + if (!isvalid(arr)) { + return str; + } + if (mxIsChar(arr)) { + char* cstr = mxArrayToString(arr); + str = std::string(cstr); + mxFree(cstr); + return str; + } + if (mxIsClass(arr,"string")) { + // convert matlab string to char and call mxArray_to_string + // TODO handle non-scalar matlab string + mxArray* arr_str = mxDuplicateArray(arr); + mxArray* arr_char; + mexCallMATLAB(1, &arr_char, 1, &arr_str, "char"); + str = mxArray_to_string(arr_char); + mxDestroyArray(arr_str); + mxDestroyArray(arr_char); + return str; + } + return str; + } + + mxArray* string_to_mxArray(const std::string& str) { + return mxCreateString(str.c_str()); + } + + MxStruct mxArray_to_mxStruct(const mxArray* s) { + if (!s) return MxStruct(); + MxStruct result; + const int nFields = mxGetNumberOfFields(s); + for (int i = 0; i < nFields; ++i) { + const char* fieldName = mxGetFieldNameByNumber(s, i); + mxArray* fieldValue = mxGetFieldByNumber(s, 0, i); + result.insert(fieldName, fieldValue); + } + return result; + } + + mxArray* mxStruct_to_mxArray(const MxStruct& mxStruct) { + Vector fieldNames; + for (const auto& [key, value]: mxStruct) { + fieldNames.push_back(key.c_str()); + } + mxArray* s = mxCreateStructMatrix(1, 1, static_cast(fieldNames.size()), fieldNames.data()); + for (const auto& [key, value]: mxStruct) { + mxSetField(s, 0, key.c_str(), value); + } + return s; + } + + MxStruct options_to_mxStruct(const Options& uno_options) { + MxStruct mxStruct; + // iterate over options + for (const auto& [option_name, option_type]: Options::option_types) { + try { + switch (option_type) { + case OptionType::INTEGER: { + const double option_value = static_cast(uno_options.get_int(option_name)); + mxStruct.insert(option_name, scalar_to_mxArray(option_value)); + break; + } + case OptionType::DOUBLE: { + const double option_value = uno_options.get_double(option_name); + mxStruct.insert(option_name, scalar_to_mxArray(option_value)); + break; + } + case OptionType::BOOL: { + const bool option_value = uno_options.get_bool(option_name); + mxStruct.insert(option_name, scalar_to_mxArray(option_value)); + break; + } + case OptionType::STRING: { + const std::string& option_value = uno_options.get_string(option_name); + mxStruct.insert(option_name, string_to_mxArray(option_value)); + break; + } + } + } catch (const std::out_of_range&) { + continue; + } + } + return mxStruct; + } + + void mxStruct_to_options(const MxStruct& options, Options& uno_options) { + for (const auto& [field_name, field_value]: options) { + try { + // set option with check between type and mxArray class + const OptionType option_type = uno_options.get_option_type(field_name); + // OptionType::DOUBLE accepts double + if (option_type == OptionType::DOUBLE) { + if (isa(field_value)) { + uno_options.set_double(field_name, mxArray_to_scalar(field_value)); + } + else { + INFO << ErrorString::format_error(ErrorType::OPTION_TYPE, field_name.c_str()) << std::endl; + } + } + // OptionType::INTEGER accepts double + else if (option_type == OptionType::INTEGER) { + if (isa(field_value)) { + uno_options.set_integer(field_name, static_cast(mxArray_to_scalar(field_value))); // cast to int + } + else { + INFO << ErrorString::format_error(ErrorType::OPTION_TYPE, field_name.c_str()) << std::endl; + } + } + // OptionType::BOOL accepts logical + else if (option_type == OptionType::BOOL) { + if (isa(field_value)) { + uno_options.set_bool(field_name, mxArray_to_scalar(field_value)); + } + else { + INFO << ErrorString::format_error(ErrorType::OPTION_TYPE, field_name.c_str()) << std::endl; + } + } + // OptionType::STRING accepts char + else if (option_type == OptionType::STRING) { + if (isa(field_value) || isa(field_value)) { + uno_options.set_string(field_name, mxArray_to_string(field_value)); + } + else { + INFO << ErrorString::format_error(ErrorType::OPTION_TYPE, field_name.c_str()) << std::endl; + } + } + } catch (const std::out_of_range&) { + // set the option with type depending on mxArray class (only double, bool, char) + if (isa(field_value)) { + uno_options.set_double(field_name, mxArray_to_scalar(field_value)); + } + else if (isa(field_value)) { + uno_options.set_bool(field_name, mxArray_to_scalar(field_value)); + } + else if (isa(field_value) || isa(field_value)) { + uno_options.set_string(field_name, mxArray_to_string(field_value)); + } + else { + INFO << ErrorString::format_error(ErrorType::OPTION_TYPE, field_name.c_str()) << std::endl; + } + } + } + } + + MxStruct result_to_mxStruct(const Result& uno_result) { + MxStruct result; + result.insert("optimization_status", scalar_to_mxArray(static_cast(uno_result.optimization_status))); + result.insert("solution_status", scalar_to_mxArray(static_cast(uno_result.solution_status))); + result.insert("solution_objective", scalar_to_mxArray(uno_result.solution_objective)); + result.insert("solution_primal_feasibility", scalar_to_mxArray(uno_result.solution_primal_feasibility)); + result.insert("solution_stationarity", scalar_to_mxArray(uno_result.solution_stationarity)); + result.insert("solution_complementarity", scalar_to_mxArray(uno_result.solution_complementarity)); + result.insert("primal_solution", vector_to_mxArray(uno_result.primal_solution)); + result.insert("constraint_dual_solution", vector_to_mxArray(uno_result.constraint_dual_solution)); + result.insert("lower_bound_dual_solution", vector_to_mxArray(uno_result.lower_bound_dual_solution)); + result.insert("upper_bound_dual_solution", vector_to_mxArray(uno_result.upper_bound_dual_solution)); + result.insert("constraint_values", vector_to_mxArray(uno_result.constraint_values)); + result.insert("number_iterations", scalar_to_mxArray(static_cast(uno_result.number_iterations))); + result.insert("number_variables", scalar_to_mxArray(static_cast(uno_result.number_variables))); + result.insert("number_constraints", scalar_to_mxArray(static_cast(uno_result.number_constraints))); + result.insert("base_indexing", scalar_to_mxArray(static_cast(uno_result.base_indexing))); + result.insert("cpu_time", scalar_to_mxArray(uno_result.cpu_time)); + result.insert("number_objective_evaluations", scalar_to_mxArray(static_cast(uno_result.number_objective_evaluations))); + result.insert("number_constraint_evaluations", scalar_to_mxArray(static_cast(uno_result.number_constraint_evaluations))); + result.insert("number_objective_gradient_evaluations", scalar_to_mxArray(static_cast(uno_result.number_objective_gradient_evaluations))); + result.insert("number_jacobian_evaluations", scalar_to_mxArray(static_cast(uno_result.number_jacobian_evaluations))); + result.insert("number_hessian_evaluations", scalar_to_mxArray(static_cast(uno_result.number_hessian_evaluations))); + result.insert("number_subproblems_solved", scalar_to_mxArray(static_cast(uno_result.number_subproblems_solved))); + return result; + } + +} // namespace \ No newline at end of file diff --git a/interfaces/Matlab/unomex/unomex_conversion.hpp b/interfaces/Matlab/unomex/unomex_conversion.hpp new file mode 100644 index 000000000..b78f6e96b --- /dev/null +++ b/interfaces/Matlab/unomex/unomex_conversion.hpp @@ -0,0 +1,99 @@ +// Copyright (c) 2025 Stefano Lovato and Charlie Vanaret +// Licensed under the MIT license. See LICENSE file in the project directory for details. + +#ifndef UNOMEX_CONVERSION_H +#define UNOMEX_CONVERSION_H + +#include +#include "linear_algebra/Vector.hpp" +#include "options/Options.hpp" +#include "optimization/Result.hpp" +#include "../cpp_classes/MxStruct.hpp" +#include "unomex_utils.hpp" +#include "mex.h" + +namespace uno { + + template + void mxArray_to_pointer(const mxArray* arr, T* pointer) { + const size_t n = mxGetNumberOfElements(arr); + const T* ptr = static_cast(mxGetData(arr)); + std::copy(ptr, ptr+n, pointer); + } + + template + mxArray* pointer_to_mxArray(const T* pointer, const size_t n) { + const mxClassID classid = get_mxClassID(); + mxArray* arr = mxCreateNumericMatrix(n, 1, classid, mxREAL); + T* ptr = static_cast(mxGetData(arr)); + std::copy(pointer, pointer+n, ptr); + return arr; + } + + template + T mxArray_to_scalar(const mxArray* arr) { + return *(static_cast(mxGetData(arr))); + } + + template + mxArray* scalar_to_mxArray(const T value) { + const mxClassID classid = get_mxClassID(); + mxArray* arr = mxCreateNumericMatrix(1, 1, classid, mxREAL); + T* ptr = static_cast(mxGetData(arr)); + *ptr = value; + return arr; + } + + template + void mxArray_to_vector(const mxArray* arr, Vector& vec) { + const size_t n = mxGetNumberOfElements(arr); // assume vec has the correct size + const T* ptr = mxGetPr(arr); + std::copy(ptr, ptr+n, vec.data()); + } + template + Vector mxArray_to_vector(const mxArray* arr) { + const size_t n = mxGetNumberOfElements(arr); + Vector vec(n); + mxArray_to_vector(arr, vec); + return vec; + } + + template + std::vector mxArray_to_stdvector(const mxArray* arr) { + const size_t n = mxGetNumberOfElements(arr); + const T* ptr = mxGetPr(arr); + std::vector vec(n); + std::copy(ptr, ptr+n, vec.data()); + return vec; + } + + template + mxArray* vector_to_mxArray(const Vector& vec) { + const mwSize n = static_cast(vec.size()); + const mxClassID classid = get_mxClassID(); + mxArray* arr = mxCreateNumericMatrix(n, 1, classid, mxREAL); + T* ptr = static_cast(mxGetData(arr)); + std::copy(vec.begin(), vec.end(), ptr); + return arr; + } + template + mxArray* vector_to_mxArray(const Vector& vec, const size_t n) { + const mxClassID classid = get_mxClassID(); + mxArray* arr = mxCreateNumericMatrix(static_cast(n), 1, classid, mxREAL); + T* ptr = static_cast(mxGetData(arr)); + std::copy(vec.begin(), vec.begin()+n, ptr); + return arr; + } + + std::string mxArray_to_string(const mxArray* arr); + mxArray* string_to_mxArray(const std::string& str); + + MxStruct mxArray_to_mxStruct(const mxArray* s); + mxArray* mxStruct_to_mxArray(const MxStruct& mxStruct); + MxStruct options_to_mxStruct(const Options& uno_options); + void mxStruct_to_options(const MxStruct& options, Options& uno_options); + MxStruct result_to_mxStruct(const Result& uno_result); + +} // namespace + +#endif // UNOMEX_CONVERSION_H \ No newline at end of file diff --git a/interfaces/Matlab/unomex/unomex_function.cpp b/interfaces/Matlab/unomex/unomex_function.cpp new file mode 100644 index 000000000..0b3bd5367 --- /dev/null +++ b/interfaces/Matlab/unomex/unomex_function.cpp @@ -0,0 +1,50 @@ +// Copyright (c) 2025 Stefano Lovato and Charlie Vanaret +// Licensed under the MIT license. See LICENSE file in the project directory for details. + +#include "unomex_function.hpp" +#include "unomex_conversion.hpp" +#include "../cpp_classes/ErrorString.hpp" + +namespace uno { + + MatlabFunctionError::MatlabFunctionError(const std::string& str) + : msg(ErrorString::format_error(ErrorType::EVALUATION) + "\n" + str) {} + + const char* MatlabFunctionError::what() const noexcept { + return this->msg.c_str(); + } + + void call_matlab_function(handle_t handle, const std::vector& inputs, std::vector& outputs) { + std::vector inputs_all = inputs; + inputs_all.insert(inputs_all.begin(), handle); + const int nrhs = static_cast(inputs_all.size()); + const int nlhs = static_cast(outputs.size()); + mxArray** prhs = inputs_all.data(); + mxArray** plhs = outputs.data(); + mxArray* err = mexCallMATLABWithTrap(nlhs, plhs, nrhs, prhs, "feval"); + if (err) { + mxArray* errmsg; + mexCallMATLAB(1, &errmsg, 1, &err, "getReport"); + const std::string strerr = mxArrayToString(errmsg); + mxDestroyArray(err); + mxDestroyArray(errmsg); + throw MatlabFunctionError(strerr); + } + } + + int32_t nargin(handle_t handle) { + mxArray* out; + mexCallMATLAB(1, &out, 1, &handle, "nargin"); + const double n = mxArray_to_scalar(out); + mxDestroyArray(out); + return static_cast(n); + } + + int32_t nargout(handle_t handle) { + mxArray* out; + mexCallMATLAB(1, &out, 1, &handle, "nargout"); + const double n = mxArray_to_scalar(out); + mxDestroyArray(out); + return static_cast(n); + } +} \ No newline at end of file diff --git a/interfaces/Matlab/unomex/unomex_function.hpp b/interfaces/Matlab/unomex/unomex_function.hpp new file mode 100644 index 000000000..79d70ea11 --- /dev/null +++ b/interfaces/Matlab/unomex/unomex_function.hpp @@ -0,0 +1,28 @@ +// Copyright (c) 2025 Stefano Lovato and Charlie Vanaret +// Licensed under the MIT license. See LICENSE file in the project directory for details. + +#ifndef UNOMEX_FUNCTION_H +#define UNOMEX_FUNCTION_H + +#include +#include +#include "optimization/EvaluationErrors.hpp" +#include "mex.h" + +namespace uno { + + // matlab function handle type + typedef mxArray* handle_t; + + struct MatlabFunctionError final : EvaluationError { + std::string msg; + explicit MatlabFunctionError(const std::string& str); + [[nodiscard]] const char* what() const noexcept override; + }; + + void call_matlab_function(handle_t handle, const std::vector& inputs, std::vector& outputs); + +} // namespace + + +#endif // UNOMEX_FUNCTION_H \ No newline at end of file diff --git a/interfaces/Matlab/unomex/unomex_utils.cpp b/interfaces/Matlab/unomex/unomex_utils.cpp new file mode 100644 index 000000000..19066af11 --- /dev/null +++ b/interfaces/Matlab/unomex/unomex_utils.cpp @@ -0,0 +1,80 @@ +// Copyright (c) 2025 Stefano Lovato and Charlie Vanaret +// Licensed under the MIT license. See LICENSE file in the project directory for details. + +#include +#include +#include "unomex_utils.hpp" +#include "unomex_function.hpp" +#include "../cpp_classes/MxStruct.hpp" + +namespace uno { + + MxArrayVectorGuard::MxArrayVectorGuard(const std::vector& vector) + : vector(vector) {} + + MxArrayVectorGuard::~MxArrayVectorGuard() { + for (auto* arr : vector) { + mxDestroyArray(arr); + } + } + + template<> mxClassID get_mxClassID() { return mxDOUBLE_CLASS; } + template<> mxClassID get_mxClassID() { return mxSINGLE_CLASS; } + template<> mxClassID get_mxClassID() { return mxINT8_CLASS; } + template<> mxClassID get_mxClassID() { return mxUINT8_CLASS; } + template<> mxClassID get_mxClassID() { return mxINT16_CLASS; } + template<> mxClassID get_mxClassID() { return mxUINT16_CLASS; } + template<> mxClassID get_mxClassID() { return mxINT32_CLASS; } + template<> mxClassID get_mxClassID() { return mxUINT32_CLASS; } + template<> mxClassID get_mxClassID() { return mxINT64_CLASS; } + template<> mxClassID get_mxClassID() { return mxUINT64_CLASS; } + template<> mxClassID get_mxClassID() { return mxLOGICAL_CLASS; } + template<> mxClassID get_mxClassID() { return mxCHAR_CLASS; } + template<> mxClassID get_mxClassID() { return mxVOID_CLASS; } + template<> mxClassID get_mxClassID() { return mxFUNCTION_CLASS; } + template<> mxClassID get_mxClassID() { return mxSTRUCT_CLASS; } + template<> mxClassID get_mxClassID() { return static_cast(mxSTRING_CLASS); } + + bool isvalid(const mxArray* arr) { + return (arr != nullptr); + } + + bool isempty(const mxArray* arr) { + return mxIsEmpty(arr); + } + + bool has_size(const mxArray* arr, const size_t nrows, const size_t ncolumns) { + if (mxGetNumberOfDimensions(arr)>2) { + return false;; + } + return (static_cast(mxGetM(arr))==nrows) && (static_cast(mxGetN(arr))==ncolumns); + } + + bool ispositive(const mxArray* arr) { + const double value = mxGetScalar(arr); + return value>=0; + } + + bool isinteger(const mxArray* arr) { + const double value = mxGetScalar(arr); + if (!std::isfinite(value)) { + return false; + } + double intpart; + return std::modf(value, &intpart) == 0.0; + } + + bool isunitary(const mxArray* arr) { + const double value = mxGetScalar(arr); + return (value==1.0) || (value==-1.0); + } + + bool isscalar(const mxArray* arr) { + return mxGetNumberOfElements(arr)==1; + } + + bool isvector(const mxArray* arr, const size_t len) { + return has_size(arr, len, 1) || has_size(arr, 1, len); + } + +}; // namespace \ No newline at end of file diff --git a/interfaces/Matlab/unomex/unomex_utils.hpp b/interfaces/Matlab/unomex/unomex_utils.hpp new file mode 100644 index 000000000..723b16816 --- /dev/null +++ b/interfaces/Matlab/unomex/unomex_utils.hpp @@ -0,0 +1,52 @@ +// Copyright (c) 2025 Stefano Lovato and Charlie Vanaret +// Licensed under the MIT license. See LICENSE file in the project directory for details. + +#ifndef UNOMEX_UTILS_H +#define UNOMEX_UTILS_H + +#include +#include +#include +#include "linear_algebra/Vector.hpp" +#include "mex.h" + +namespace uno { + + // guard on std::vector for automatic memory management + struct MxArrayVectorGuard { + const std::vector& vector; + explicit MxArrayVectorGuard(const std::vector& vector); + ~MxArrayVectorGuard(); + }; + + constexpr int mxSTRING_CLASS = 19; // cf. https://it.mathworks.com/matlabcentral/answers/2102376-how-do-i-process-a-string-class-in-a-mex-function + + template + mxClassID get_mxClassID(); + + bool isvalid(const mxArray* arr); + bool isempty(const mxArray* arr); + + template + bool isa(const mxArray* arr) { + return mxGetClassID(arr) == get_mxClassID(); + } + + bool has_size(const mxArray* arr, size_t nrows, size_t ncolumns); + bool ispositive(const mxArray* arr); + bool isinteger(const mxArray* arr); + bool isunitary(const mxArray* arr); + bool isscalar(const mxArray* arr); + bool isvector(const mxArray* arr, size_t len); + + template + Vector convert_vector_type(const Vector& input) { + Vector output(input.size()); + std::transform(input.begin(), input.end(), output.begin(), + [](const InType& val) { return static_cast(val); }); + return output; + } + +}; // namespace + +#endif // UNOMEX_UTILS_H \ No newline at end of file diff --git a/interfaces/Matlab/unomex/unomex_validation.cpp b/interfaces/Matlab/unomex/unomex_validation.cpp new file mode 100644 index 000000000..7af43857f --- /dev/null +++ b/interfaces/Matlab/unomex/unomex_validation.cpp @@ -0,0 +1,169 @@ +// Copyright (c) 2025 Stefano Lovato and Charlie Vanaret +// Licensed under the MIT license. See LICENSE file in the project directory for details. + +#include "../cpp_classes/ErrorString.hpp" +#include "unomex_utils.hpp" +#include "unomex_conversion.hpp" +#include "unomex_validation.hpp" + +namespace uno { + + // validate Matlab function handle + bool validate_matlab_handle_field(handle_t handle, const std::string& field_name, std::string& errmsg) { + // Matlab throws error when wrong nargin/nargour etc... just check the validity of the handle + if (!isvalid(handle)) { + errmsg = ErrorString::format_error(ErrorType::MISSING_FIELD, field_name.c_str()); + return false; + } + if (!isa(handle)) { + errmsg = ErrorString::format_error(ErrorType::INVALID_HANDLE); + return false; + } + return true; + } + + bool validate_struct_input(const mxArray* arr, const size_t position, std::string& errmsg) { + if (!isvalid(arr) || !isa(arr)) { + errmsg = ErrorString::format_error(ErrorType::INPUT_STRUCT, position); + return false; + } + if (!isscalar(arr)) { + errmsg = ErrorString::format_error(ErrorType::INPUT_SCALAR, position); + return false; + } + return true; + } + + bool validate_char_field(const mxArray* arr, const std::vector& chars, const std::string& field_name, std::string& errmsg) { + if (!isvalid(arr)) { + errmsg = ErrorString::format_error(ErrorType::MISSING_FIELD, field_name.c_str()); + return false; + } + if (!isa(arr) || !isscalar(arr)) { + errmsg = ErrorString::format_error(ErrorType::FIELD_CHAR, field_name.c_str()); + return false; + } + if (!chars.empty()) { + if (const char c = mxArray_to_scalar(arr); + std::find(chars.begin(), chars.end(), c) == chars.end()) { + errmsg = "Field '" + field_name + "' must be one of { "; + for (const size_t i: Range(chars.size())) { + if (i > 0) { + errmsg += ", "; + } + errmsg += "'" + std::string({chars[i]}) + "'"; + } + errmsg += " }."; + return false; + } + } + return true; + } + + bool validate_string_field(const mxArray* arr, const std::vector& strings, const std::string& field_name, std::string& errmsg) { + if (!isvalid(arr)) { + errmsg = ErrorString::format_error(ErrorType::MISSING_FIELD, field_name.c_str()); + return false; + } + if (!isa(arr) && !isa(arr)) { + errmsg = ErrorString::format_error(ErrorType::FIELD_STRING, field_name.c_str()); + return false; + } + if (!strings.empty()) { + if (const std::string str = mxArray_to_string(arr); + std::find(strings.begin(), strings.end(), str) == strings.end()) { + errmsg = "Field '" + field_name + "' must be one of { "; + for (const size_t i: Range(strings.size())) { + if (i > 0) { + errmsg += ", "; + } + errmsg += "'" + strings[i] + "'"; + } + errmsg += " }."; + return false; + } + } + return true; + } + + bool validate_positive_integer_field(const mxArray* arr, const std::string& field_name, std::string& errmsg) { + if (!isvalid(arr)) { + errmsg = ErrorString::format_error(ErrorType::MISSING_FIELD, field_name.c_str()); + return false; + } + if (!isa(arr) || !isscalar(arr) || !ispositive(arr) || !isinteger(arr) ) { + errmsg = ErrorString::format_error(ErrorType::FIELD_POSITIVE_INTEGER, field_name.c_str()); + return false; + } + return true; + } + + bool validate_unitary_field(const mxArray* arr, const std::string& field_name, std::string& errmsg) { + if (!isvalid(arr)) { + errmsg = ErrorString::format_error(ErrorType::MISSING_FIELD, field_name.c_str()); + return false; + } + if (!isa(arr) || !isscalar(arr) || !isinteger(arr) || !isunitary(arr) ) { + errmsg = ErrorString::format_error(ErrorType::FIELD_UNITARY, field_name.c_str()); + return false; + } + return true; + } + + bool validate_double_vector_field(const mxArray* arr, const size_t len, const std::string& field_name, std::string& errmsg) { + if (!isvalid(arr)) { + errmsg = ErrorString::format_error(ErrorType::MISSING_FIELD, field_name.c_str()); + return false; + } + if (!isa(arr) || !isvector(arr, len) ) { + errmsg = ErrorString::format_error(ErrorType::FIELD_VECTOR, field_name.c_str(), len); + return false; + } + return true; + } + + bool validate_integer_vector_field(const mxArray* arr, const size_t len, const std::string& field_name, std::string& errmsg) { + if (!isvalid(arr)) { + errmsg = ErrorString::format_error(ErrorType::MISSING_FIELD, field_name.c_str()); + return false; + } + if (!isa(arr) || !isvector(arr, len) ) { + errmsg = ErrorString::format_error(ErrorType::FIELD_VECTOR_INT, field_name.c_str(), len); + return false; + } + const Vector vector = mxArray_to_vector(arr); + double intpart; + for (auto value : vector) { + if (std::modf(value, &intpart) != 0.0) { + errmsg = ErrorString::format_error(ErrorType::FIELD_VECTOR_INT, field_name.c_str(), len); + return false; + } + } + return true; + } + + bool validate_double_scalar_output(const mxArray* arr, std::string& errmsg) { + if (!isa(arr) || !isscalar(arr)) { + errmsg = ErrorString::format_error(ErrorType::OUTPUT_SCALAR, 1); + return false; + } + return true; + } + + bool validate_double_vector_output(const mxArray* arr, size_t n, std::string& errmsg) { + if (!isa(arr) || !isvector(arr, n)) { + errmsg = ErrorString::format_error(ErrorType::OUTPUT_VECTOR, 1, n); + return false; + } + return true; + } + + bool validate_bool_scalar_output(const mxArray* arr, std::string& errmsg) { + if (!isa(arr) || !isscalar(arr)) { + errmsg = ErrorString::format_error(ErrorType::OUTPUT_BOOL_SCALAR, 1); + return false; + } + return true; + } + +} // namespace \ No newline at end of file diff --git a/interfaces/Matlab/unomex/unomex_validation.hpp b/interfaces/Matlab/unomex/unomex_validation.hpp new file mode 100644 index 000000000..1cde53778 --- /dev/null +++ b/interfaces/Matlab/unomex/unomex_validation.hpp @@ -0,0 +1,36 @@ +// Copyright (c) 2025 Stefano Lovato and Charlie Vanaret +// Licensed under the MIT license. See LICENSE file in the project directory for details. + +#ifndef UNOMEX_VALIDATION_H +#define UNOMEX_VALIDATION_H + +#include +#include +#include "unomex_function.hpp" +#include "mex.h" + +namespace uno { + + bool validate_matlab_handle_field(handle_t handle, const std::string& field_name, std::string& errmsg); + + bool validate_struct_input(const mxArray* arr, size_t position, std::string& errmsg); + + bool validate_char_field(const mxArray* arr, const std::vector& chars, const std::string& field_name, std::string& errmsg); + bool validate_string_field(const mxArray* arr, const std::vector& chars, const std::string& field_name, std::string& errmsg); + bool validate_positive_integer_field(const mxArray* arr, const std::string& field_name, std::string& errmsg); + bool validate_unitary_field(const mxArray* arr, const std::string& field_name, std::string& errmsg); + + bool validate_double_vector_field(const mxArray* arr, size_t len, const std::string& field_name, std::string& errmsg); + bool validate_integer_vector_field(const mxArray* arr, size_t len, const std::string& field_name, std::string& errmsg); + + bool validate_double_scalar_output(const mxArray* arr, std::string& errmsg); + bool validate_bool_scalar_output(const mxArray* arr, std::string& errmsg); + bool validate_double_vector_output(const mxArray* arr, size_t n, std::string& errmsg); + bool validate_gradient(const mxArray* arr, size_t n, std::string& errmsg); + bool validate_constraint_jacobian(const mxArray* arr, size_t n, std::string& errmsg); + +} // namespace + + + +#endif // UNOMEX_VALIDATION_H \ No newline at end of file diff --git a/interfaces/Python/cpp_classes/UnoSolverWrapper.hpp b/interfaces/Python/cpp_classes/UnoSolverWrapper.hpp index fa406ca4e..4de8d9196 100644 --- a/interfaces/Python/cpp_classes/UnoSolverWrapper.hpp +++ b/interfaces/Python/cpp_classes/UnoSolverWrapper.hpp @@ -4,6 +4,7 @@ #ifndef UNO_UNOSOLVERWRAPPER_H #define UNO_UNOSOLVERWRAPPER_H +#include #include #include "Uno.hpp" #include "options/Options.hpp" @@ -51,4 +52,4 @@ namespace uno { }; } // namespace -#endif // UNO_UNOSOLVERWRAPPER_H \ No newline at end of file +#endif // UNO_UNOSOLVERWRAPPER_H diff --git a/interfaces/UserModel.hpp b/interfaces/UserModel.hpp index e495d74c9..d6eb43ba3 100644 --- a/interfaces/UserModel.hpp +++ b/interfaces/UserModel.hpp @@ -4,9 +4,12 @@ #ifndef UNO_USERMODEL_H #define UNO_USERMODEL_H +#include +#include #include #include #include +#include #include "C/Uno_C_API.h" #include "optimization/ProblemType.hpp" @@ -82,6 +85,111 @@ namespace uno { DoubleVector initial_primal_iterate{}; DoubleVector initial_dual_iterate{}; }; + + // base class for user stream callback + class UserStreamCallback { + public: + UserStreamCallback() = default; + virtual ~UserStreamCallback() = default; + virtual int32_t operator()(const char* buf, int32_t len) const = 0; + }; + + // std::streambuf wrapper around UserStreamCallback + class UserStreamBuffer : public std::streambuf { + public: + UserStreamBuffer(UserStreamCallback* user_stream_callback, std::size_t buffer_size) : + user_stream_callback(user_stream_callback) { + // allocate output buffer and set stream buffer pointer + this->buffer = new char[buffer_size]; + this->setp(this->buffer, this->buffer + buffer_size - 1); + } + ~UserStreamBuffer() override { + // flush remaining data and release buffer memory + this->sync(); + // user_stream_callback has to be deleted explictly + delete[] this->buffer; + } + + protected: + // called on buffer overflow + int overflow(int character = EOF) override { + if (character != EOF) { + // insert the character into the buffer + *this->pptr() = traits_type::to_char_type(character); + this->pbump(1); + } + // return EOF for error + return (this->flush_buffer() == 0) ? character : EOF; + } + + // put char into buffer with \n triggering flush + std::streamsize xsputn(const char* s, std::streamsize n) override { + const char* p = s; + const char* end = s + n; + while (p < end) { + const char* pos = static_cast(std::memchr(p, '\n', end - p)); + if (pos) { + std::ptrdiff_t chunk_len = pos - p + 1; + if (this->epptr() - this->pptr() < chunk_len) + this->overflow(); + std::memcpy(this->pptr(), p, chunk_len); + this->pbump(static_cast(chunk_len)); + this->overflow(); + p = pos + 1; + } else { + std::ptrdiff_t chunk_len = end - p; + if (this->epptr() - this->pptr() < chunk_len) { + overflow(); + if (chunk_len > (this->epptr() - this->pptr())) + chunk_len = this->epptr() - this->pptr(); + } + std::memcpy(pptr(), p, chunk_len); + this->pbump(static_cast(chunk_len)); + p += chunk_len; + } + } + return n; + } + + int sync() override { + return this->flush_buffer(); + } + + private: + UserStreamCallback* user_stream_callback; + char* buffer; + + // flush buffer to the user_stream_callback + int flush_buffer() { + // check for invalid stream callback + if (!this->user_stream_callback) { + return -1; + } + std::ptrdiff_t current_used_buffer_size = this->pptr() - this->pbase(); + if (current_used_buffer_size > 0) { + // call user stream callback + const int32_t callback_result = (*this->user_stream_callback)(this->pbase(), static_cast(current_used_buffer_size)); + if (callback_result != static_cast(current_used_buffer_size)) { + return -1; + } + // move buffer pointer + this->pbump(static_cast(-current_used_buffer_size)); + } + return 0; + } + }; + + // std::ostream wrapper around UserStreamCallback + class UserOStream : public std::ostream { + public: + UserOStream(UserStreamCallback* user_stream_callback, std::size_t buffer_size = 1024) : // 1024 default buffer size, sufficient to store the entire line + std::ostream(&this->buffer), buffer(user_stream_callback, buffer_size) { } + + private: + // internal stream buffer that sends output to the LoggerStreamUserCallback + UserStreamBuffer buffer; + }; + } // namespace #endif // UNO_USERMODEL_H \ No newline at end of file diff --git a/pyproject.toml b/pyproject.toml index 3c553be61..375956d19 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -35,7 +35,7 @@ wheel.packages = ["interfaces/Python/unopy"] # point to a dir containing __init # support for passing paths to dependencies from environment to CMake [tool.scikit-build.cmake.define] -BLAS_LIBRARIES = "dependencies/lib/libblas.a;dependencies/lib/libcblas.a" +BLAS_LIBRARIES = "dependencies/lib/libcblas.a;dependencies/lib/libblas.a" LAPACK_LIBRARIES = "dependencies/lib/liblapack.a" # BQPD (cmake -DBQPD=/path/.../libbqpd.a) BQPD = "dependencies/lib/libbqpd.a" diff --git a/uno/options/Options.cpp b/uno/options/Options.cpp index 4f0b8462b..41caa5b91 100644 --- a/uno/options/Options.cpp +++ b/uno/options/Options.cpp @@ -339,4 +339,4 @@ namespace uno { DISCRETE << "\nNon-default options:\n" << option_list << '\n'; } } -} // namespace +} // namespace \ No newline at end of file