From 9415957080e645d8eafa9c05cc7926c1e5ad314b Mon Sep 17 00:00:00 2001 From: vcantarella Date: Mon, 23 Jun 2025 15:43:20 +0200 Subject: [PATCH 1/2] added honour_x0 interface for the algorithm bobyqa in the C api --- c/bobyqa_c.f90 | 13 +++++++----- c/examples/bobyqa/bobyqa_example.c | 32 +++++++++++++++++++++++++++++- c/include/prima/prima.h | 4 ++++ c/prima.c | 8 ++++++-- 4 files changed, 49 insertions(+), 8 deletions(-) diff --git a/c/bobyqa_c.f90 b/c/bobyqa_c.f90 index 6f318ec742..dc810c4fea 100644 --- a/c/bobyqa_c.f90 +++ b/c/bobyqa_c.f90 @@ -14,8 +14,8 @@ module bobyqa_c_mod subroutine bobyqa_c(cobj_ptr, data_ptr, n, x, f, xl, xu, nf, rhobeg, rhoend, & - & ftarget, maxfun, npt, iprint, callback_ptr, info) bind(C) -use, intrinsic :: iso_c_binding, only : C_DOUBLE, C_INT, C_FUNPTR, C_PTR, C_ASSOCIATED, C_F_PROCPOINTER, C_F_POINTER + & ftarget, maxfun, npt, iprint, honour_x0, callback_ptr, info) bind(C) +use, intrinsic :: iso_c_binding, only : C_DOUBLE, C_INT, C_FUNPTR, C_PTR, C_ASSOCIATED, C_F_PROCPOINTER, C_F_POINTER, C_BOOL use, non_intrinsic :: bobyqa_mod, only : bobyqa use, non_intrinsic :: cintrf_mod, only : COBJ, CCALLBACK use, non_intrinsic :: consts_mod, only : RP, IK @@ -38,6 +38,7 @@ subroutine bobyqa_c(cobj_ptr, data_ptr, n, x, f, xl, xu, nf, rhobeg, rhoend, & integer(C_INT), intent(in), value :: maxfun integer(C_INT), intent(in), value :: npt integer(C_INT), intent(in), value :: iprint +logical(C_BOOL), intent(in), value :: honour_x0 type(C_FUNPTR), intent(in), value :: callback_ptr integer(C_INT), intent(out) :: info @@ -62,7 +63,7 @@ subroutine bobyqa_c(cobj_ptr, data_ptr, n, x, f, xl, xu, nf, rhobeg, rhoend, & real(RP), allocatable :: rhoend_loc real(RP), allocatable :: xl_loc(:) real(RP), allocatable :: xu_loc(:) - +logical:: honour_x0_loc ! Read the inputs and convert them to the Fortran side types ! The following inputs correspond to compulsory arguments in the Fortran code. @@ -100,6 +101,8 @@ subroutine bobyqa_c(cobj_ptr, data_ptr, n, x, f, xl, xu, nf, rhobeg, rhoend, & npt_loc = int(npt, kind(npt_loc)) end if iprint_loc = int(iprint, kind(iprint_loc)) +honour_x0_loc = logical(honour_x0, kind(.true.)) + ! Call the Fortran code if (c_associated(callback_ptr)) then @@ -108,10 +111,10 @@ subroutine bobyqa_c(cobj_ptr, data_ptr, n, x, f, xl, xu, nf, rhobeg, rhoend, & call c_f_procpointer(callback_ptr, cb_ptr) ! We then provide the closure to the algorithm. call bobyqa(calfun, x_loc, f_loc, xl=xl_loc, xu=xu_loc, nf=nf_loc, rhobeg=rhobeg_loc, rhoend=rhoend_loc, & - & ftarget=ftarget_loc, maxfun=maxfun_loc, npt=npt_loc, iprint=iprint_loc, callback_fcn=callback_fcn, info=info_loc) + & ftarget=ftarget_loc, maxfun=maxfun_loc, npt=npt_loc, iprint=iprint_loc, honour_x0=honour_x0_loc, callback_fcn=callback_fcn, info=info_loc) else call bobyqa(calfun, x_loc, f_loc, xl=xl_loc, xu=xu_loc, nf=nf_loc, rhobeg=rhobeg_loc, rhoend=rhoend_loc, & - & ftarget=ftarget_loc, maxfun=maxfun_loc, npt=npt_loc, iprint=iprint_loc, info=info_loc) + & ftarget=ftarget_loc, maxfun=maxfun_loc, npt=npt_loc, iprint=iprint_loc, honour_x0=honour_x0_loc, info=info_loc) end if ! Write the outputs diff --git a/c/examples/bobyqa/bobyqa_example.c b/c/examples/bobyqa/bobyqa_example.c index b73b1dc4df..8fa4f34e19 100644 --- a/c/examples/bobyqa/bobyqa_example.c +++ b/c/examples/bobyqa/bobyqa_example.c @@ -68,5 +68,35 @@ int main(int argc, char * argv[]) // Free the result prima_free_result(&result); - return success; + // for the second part of the example, we will change the initial point close to the upper bound + // and set honour_x0 to true + double new_x0[2] = {4.4, 4.4}; + + // Set up the problem + prima_problem_t honour_problem; + prima_init_problem(&honour_problem, n); + honour_problem.x0 = new_x0; + honour_problem.calfun = &fun; + // We define an upper bound that will be active in order to demonstrate the usage. + honour_problem.xl = xl; + honour_problem.xu = xu; + + prima_options_t honour_options; + prima_init_options(&honour_options); + honour_options.iprint = PRIMA_MSG_EXIT; + honour_options.rhoend = 1e-6; + honour_options.maxfun = 500*n; + honour_options.honour_x0 = true; // change this to false to see the warning of resetting x0 + honour_options.callback = &callback; + // Call the solver again with honour_x0 = true + prima_result_t honour_result; + const prima_rc_t honour_rc = prima_minimize(PRIMA_BOBYQA, honour_problem, honour_options, &honour_result); + + // Print the result + printf("x* = {%g, %g}, rc = %d, msg = '%s', evals = %d\n", honour_result.x[0], honour_result.x[1], honour_rc, honour_result.message, honour_result.nf); + + // Check the result + int honour_success = (fabs(honour_result.x[0] - 4.5) > 2e-2 || fabs(honour_result.x[1] - 4) > 2e-2); + + return honour_success; } diff --git a/c/include/prima/prima.h b/c/include/prima/prima.h index 6755765729..a3ea3c5f62 100644 --- a/c/include/prima/prima.h +++ b/c/include/prima/prima.h @@ -235,6 +235,10 @@ typedef struct { // Default: NULL, which means that no callback will be called prima_callback_t callback; + // honour_x0: whether to honour the initial point x0 + // If set to true, the solver will use the initial point x0 as the starting point even if the trust region reaches the bounds. + // If set to false, the solver will potentially ignore the initial point x0 when it is too close to the bound (relative to rhobeg). + bool honour_x0; } prima_options_t; diff --git a/c/prima.c b/c/prima.c index 175b1ca19c..3cbf9e2e9c 100644 --- a/c/prima.c +++ b/c/prima.c @@ -53,6 +53,7 @@ prima_rc_t prima_init_options(prima_options_t *const options) return PRIMA_NULL_OPTIONS; memset(options, 0, sizeof(prima_options_t)); + // by initializing with 0 we already assumes honour_x0 = false options->rhobeg = NAN; // Will be interpreted by Fortran as not present options->rhoend = NAN; // Will be interpreted by Fortran as not present options->iprint = PRIMA_MSG_NONE; @@ -198,7 +199,8 @@ int cobyla_c(const int m_nlcon, const prima_objcon_t calcfc, const void *data, c const prima_callback_t callback, int *const info); int bobyqa_c(prima_obj_t calfun, const void *data, const int n, double x[], double *const f, const double xl[], const double xu[], - int *const nf, const double rhobeg, const double rhoend, const double ftarget, const int maxfun, const int npt, const int iprint, const prima_callback_t callback, int *const info); + int *const nf, const double rhobeg, const double rhoend, const double ftarget, const int maxfun, const int npt, const int iprint, + const bool honour_x0, const prima_callback_t callback, int *const info); int newuoa_c(prima_obj_t calfun, const void *data, const int n, double x[], double *const f, int *const nf, const double rhobeg, const double rhoend, const double ftarget, const int maxfun, const int npt, const int iprint, const prima_callback_t callback, int *const info); @@ -231,7 +233,9 @@ prima_rc_t prima_minimize(const prima_algorithm_t algorithm, const prima_problem switch (algorithm) { case PRIMA_BOBYQA: - bobyqa_c(problem.calfun, options.data, problem.n, result->x, &(result->f), problem.xl, problem.xu, &(result->nf), options.rhobeg, options.rhoend, options.ftarget, options.maxfun, options.npt, options.iprint, options.callback, &info); + bobyqa_c(problem.calfun, options.data, problem.n, result->x, &(result->f), problem.xl, problem.xu, + &(result->nf), options.rhobeg, options.rhoend, options.ftarget, options.maxfun, options.npt, options.iprint, + options.honour_x0, options.callback, &info); result->cstrv = 0.0; break; From f47c3d813b21442f01c33bf1f75ab45bffe4d4fb Mon Sep 17 00:00:00 2001 From: Grant Hecht Date: Sun, 13 Jul 2025 11:49:04 -0400 Subject: [PATCH 2/2] Continue long argument line to new line to avoid line truncated error --- c/bobyqa_c.f90 | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/c/bobyqa_c.f90 b/c/bobyqa_c.f90 index dc810c4fea..d9fb9f4130 100644 --- a/c/bobyqa_c.f90 +++ b/c/bobyqa_c.f90 @@ -111,10 +111,12 @@ subroutine bobyqa_c(cobj_ptr, data_ptr, n, x, f, xl, xu, nf, rhobeg, rhoend, & call c_f_procpointer(callback_ptr, cb_ptr) ! We then provide the closure to the algorithm. call bobyqa(calfun, x_loc, f_loc, xl=xl_loc, xu=xu_loc, nf=nf_loc, rhobeg=rhobeg_loc, rhoend=rhoend_loc, & - & ftarget=ftarget_loc, maxfun=maxfun_loc, npt=npt_loc, iprint=iprint_loc, honour_x0=honour_x0_loc, callback_fcn=callback_fcn, info=info_loc) + & ftarget=ftarget_loc, maxfun=maxfun_loc, npt=npt_loc, iprint=iprint_loc, & + & honour_x0=honour_x0_loc, callback_fcn=callback_fcn, info=info_loc) else call bobyqa(calfun, x_loc, f_loc, xl=xl_loc, xu=xu_loc, nf=nf_loc, rhobeg=rhobeg_loc, rhoend=rhoend_loc, & - & ftarget=ftarget_loc, maxfun=maxfun_loc, npt=npt_loc, iprint=iprint_loc, honour_x0=honour_x0_loc, info=info_loc) + & ftarget=ftarget_loc, maxfun=maxfun_loc, npt=npt_loc, iprint=iprint_loc, & + & honour_x0=honour_x0_loc, info=info_loc) end if ! Write the outputs