added honour_x0 interface for the algorithm bobyqa in the C api

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

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;
 }

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;
 
 

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;