Manually applied jschueller's closed callback PR

.github/actions/spelling/allow.txt

@@ -45,6 +45,7 @@ BQPGASIM
 BROWNAL
 Bindel
 Broyden
+CCALLBACK
 CGRAD
 CHANDHEQ
 CHCKTST
@@ -2048,6 +2049,7 @@ COBJCON
 cobjfun
 cobjfuncon
 constrc
+evalcallback
 evalcobj
 evalcobjcon
 execstack
@@ -2133,3 +2135,5 @@ TWOBARS
 lang
 archnorma
 orthtol
+fcb
+fcn

c/bobyqa_c.f90

@@ -13,7 +13,8 @@ module bobyqa_c_mod
 contains
 
 
-subroutine bobyqa_c(cobj_ptr, data_ptr, n, x, f, xl, xu, nf, rhobeg, rhoend, ftarget, maxfun, npt, iprint, info) bind(C)
+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
 use, non_intrinsic :: cintrf_mod, only : COBJ
 use, non_intrinsic :: consts_mod, only : RP, IK
@@ -36,6 +37,7 @@ subroutine bobyqa_c(cobj_ptr, data_ptr, n, x, f, xl, xu, nf, rhobeg, rhoend, fta
 integer(C_INT), intent(in), value :: maxfun
 integer(C_INT), intent(in), value :: npt
 integer(C_INT), intent(in), value :: iprint
+type(C_FUNPTR), intent(in), value :: callback_ptr
 integer(C_INT), intent(out) :: info
 
 ! Local variables
@@ -65,7 +67,7 @@ subroutine bobyqa_c(cobj_ptr, data_ptr, n, x, f, xl, xu, nf, rhobeg, rhoend, fta
 
 ! Call the Fortran code
 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, callback_fcn=callback_fcn, info=info_loc)
 
 ! Write the outputs
 x = real(x_loc, kind(x))
@@ -89,6 +91,26 @@ subroutine calfun(x_sub, f_sub)
 call evalcobj(cobj_ptr, data_ptr, x_sub, f_sub)
 end subroutine calfun
 
+
+subroutine callback_fcn(x_sub, f_sub, nf_sub, tr_sub, cstrv_sub, nlconstr_sub, terminate_sub)
+use, non_intrinsic :: consts_mod, only : RP, IK
+use, non_intrinsic :: cintrf_mod, only : evalcallback
+use, intrinsic :: iso_c_binding, only : C_DOUBLE, C_INT, C_ASSOCIATED
+implicit none
+real(RP), intent(in) :: x_sub(:)
+real(RP), intent(in) :: f_sub
+integer(IK), intent(in) :: nf_sub
+integer(IK), intent(in) :: tr_sub
+real(RP), intent(in) :: cstrv_sub
+real(RP), intent(in) :: nlconstr_sub(:)
+logical, intent(out) :: terminate_sub
+terminate_sub = .false.
+if (C_ASSOCIATED(callback_ptr)) then
+    call evalcallback(callback_ptr, x_sub, f_sub, nf_sub, tr_sub, cstrv_sub, nlconstr_sub, terminate_sub)
+end if
+end subroutine callback_fcn
+
+
 end subroutine bobyqa_c
 
 

c/cintrf.f90

@@ -7,7 +7,7 @@ module cintrf_mod
 
 implicit none
 private
-public :: COBJ, COBJCON, evalcobj, evalcobjcon
+public :: COBJ, COBJCON, evalcobj, evalcobjcon, evalcallback
 
 
 abstract interface
@@ -31,6 +31,22 @@ subroutine COBJCON(x, f, constr, data_ptr) bind(c)
     type(C_PTR), intent(in), value :: data_ptr
     end subroutine COBJCON
 
+    subroutine CCALLBACK(n, x, f, nf, tr, cstrv, m_nlcon, nlconstr, terminate) bind(c)
+    use, intrinsic :: iso_c_binding, only : C_DOUBLE, C_BOOL, C_INT
+    implicit none
+    integer(C_INT), intent(in), value :: n
+    ! We cannot use assumed-shape arrays for C interoperability
+    real(C_DOUBLE), intent(in) :: x(*)
+    real(C_DOUBLE), intent(in), value :: f
+    integer(C_INT), intent(in), value :: nf
+    integer(C_INT), intent(in), value :: tr
+    real(C_DOUBLE), intent(in), value :: cstrv
+    integer(C_INT), intent(in), value :: m_nlcon
+    real(C_DOUBLE), intent(in) :: nlconstr(*)
+    logical(C_BOOL), intent(out) :: terminate
+    end subroutine CCALLBACK
+
+
 end interface
 
 
@@ -93,4 +109,48 @@ subroutine evalcobjcon(cobjcon_ptr, data_ptr, x, f, constr)
 
 end subroutine evalcobjcon
 
+
+subroutine evalcallback(fcb_ptr, x, f, nf, tr, cstrv, nlconstr, terminate)
+use, non_intrinsic :: consts_mod, only : RP, IK
+use, intrinsic :: iso_c_binding, only : C_DOUBLE, C_FUNPTR, C_F_PROCPOINTER, C_BOOL, C_INT
+implicit none
+type(C_FUNPTR), intent(in) :: fcb_ptr
+real(RP), intent(in) :: x(:)
+real(RP), intent(in) :: f
+integer(IK), intent(in) :: nf
+integer(IK), intent(in) :: tr
+real(RP), intent(in) :: cstrv
+real(RP), intent(in) :: nlconstr(:)
+logical, intent(out) :: terminate
+
+! Local variables
+procedure(CCALLBACK), pointer :: cb_ptr
+integer(C_INT) :: n_loc
+real(C_DOUBLE) :: x_loc(size(x))
+real(C_DOUBLE) :: f_loc
+integer(C_INT) :: nf_loc
+integer(C_INT) :: tr_loc
+real(C_DOUBLE) :: cstrv_loc
+real(C_DOUBLE) :: nlconstr_loc(size(nlconstr))
+logical(C_BOOL) :: terminate_loc
+
+! Read the inputs and convert them to the types specified in CCALLBACK
+n_loc = size(x)
+x_loc = real(x, kind(x_loc))
+f_loc = real(f, kind(f_loc))
+nf_loc = int(nf, kind(nf_loc))
+tr_loc = int(tr, kind(tr_loc))
+cstrv_loc = real(cstrv, kind(cstrv_loc))
+nlconstr_loc = real(nlconstr, kind(nlconstr_loc))
+call C_F_PROCPOINTER(fcb_ptr, cb_ptr)
+
+! Call the C objective function
+call cb_ptr(n_loc, x_loc, f_loc, nf_loc, tr_loc, cstrv_loc, size(nlconstr_loc), nlconstr_loc, terminate_loc)
+
+! Write the output
+terminate = logical(terminate_loc, kind(terminate))
+
+end subroutine evalcallback
+
+
 end module cintrf_mod

c/cobyla_c.f90

@@ -12,8 +12,8 @@ module cobyla_c_mod
 contains
 
 
-subroutine cobyla_c(m_nlcon, cobjcon_ptr, data_ptr, n, x, f, cstrv, nlconstr, m_ineq, Aineq, bineq, m_eq, Aeq, beq, xl, xu, &
-    & f0, nlconstr0, nf, rhobeg, rhoend, ftarget, maxfun, iprint, info) bind(C)
+subroutine cobyla_c(m_nlcon, cobjcon_ptr, data_ptr, n, x, f, cstrv, nlconstr, m_ineq, Aineq, bineq, m_eq, Aeq, beq, &
+    & xl, xu, f0, nlconstr0, nf, rhobeg, rhoend, ftarget, maxfun, iprint, callback_ptr, info) bind(C)
 use, intrinsic :: iso_c_binding, only : C_DOUBLE, C_INT, C_FUNPTR, C_PTR
 use, non_intrinsic :: cintrf_mod, only : COBJCON
 use, non_intrinsic :: consts_mod, only : RP, IK
@@ -46,6 +46,7 @@ subroutine cobyla_c(m_nlcon, cobjcon_ptr, data_ptr, n, x, f, cstrv, nlconstr, m_
 real(C_DOUBLE), intent(in), value :: ftarget
 integer(C_INT), intent(in), value :: maxfun
 integer(C_INT), intent(in), value :: iprint
+type(C_FUNPTR), intent(in), value :: callback_ptr
 integer(C_INT), intent(out) :: info
 
 ! Local variables
@@ -94,7 +95,7 @@ subroutine cobyla_c(m_nlcon, cobjcon_ptr, data_ptr, n, x, f, cstrv, nlconstr, m_
     & Aineq=Aineq_loc, bineq=bineq_loc, Aeq=Aeq_loc, beq=beq_loc, &
     & xl=xl_loc, xu=xu_loc, f0=f0_loc, nlconstr0=nlconstr0_loc, nf=nf_loc, &
     & rhobeg=rhobeg_loc, rhoend=rhoend_loc, ftarget=ftarget_loc, maxfun=maxfun_loc, &
-    & iprint=iprint_loc, info=info_loc)
+    & iprint=iprint_loc, callback_fcn=callback_fcn, info=info_loc)
 
 ! Write the outputs
 x = real(x_loc, kind(x))
@@ -121,6 +122,26 @@ subroutine calcfc(x_sub, f_sub, constr_sub)
 call evalcobjcon(cobjcon_ptr, data_ptr, x_sub, f_sub, constr_sub)
 end subroutine calcfc
 
+
+subroutine callback_fcn(x_sub, f_sub, nf_sub, tr_sub, cstrv_sub, nlconstr_sub, terminate_sub)
+use, non_intrinsic :: consts_mod, only : RP, IK
+use, non_intrinsic :: cintrf_mod, only : evalcallback
+use, intrinsic :: iso_c_binding, only : C_DOUBLE, C_INT, C_ASSOCIATED
+implicit none
+real(RP), intent(in) :: x_sub(:)
+real(RP), intent(in) :: f_sub
+integer(IK), intent(in) :: nf_sub
+integer(IK), intent(in) :: tr_sub
+real(RP), intent(in) :: cstrv_sub
+real(RP), intent(in) :: nlconstr_sub(:)
+logical, intent(out) :: terminate_sub
+terminate_sub = .false.
+if (C_ASSOCIATED(callback_ptr)) then
+    call evalcallback(callback_ptr, x_sub, f_sub, nf_sub, tr_sub, cstrv_sub, nlconstr_sub, terminate_sub)
+end if
+end subroutine callback_fcn
+
+
 end subroutine cobyla_c
 
 

c/examples/bobyqa/bobyqa_example.c

@@ -12,6 +12,15 @@ static void fun(const double x[], double *f, const void *data)
   (void)data;
 }
 
+static void callback(int n, const double x[], double f, int nf, int tr, double cstrv, const int m_nlcon, const double nlconstr[], bool *terminate)
+{
+  (void)n;
+  printf("progress: x=[%g;%g] f=%g cstrv=%g nf=%d tr=%d\n", x[0], x[1], f, cstrv, nf, tr);
+  *terminate = 0;
+  (void)m_nlcon;
+  (void)nlconstr;
+}
+
 int main(int argc, char * argv[])
 {
   (void)argc;
@@ -27,6 +36,7 @@ int main(int argc, char * argv[])
   options.iprint = PRIMA_MSG_EXIT;
   options.rhoend= 1e-3;
   options.maxfun = 200*n;
+  options.callback = &callback;
   prima_result_t result;
   const int rc = prima_minimize(PRIMA_BOBYQA, &problem, &options, &result);
   printf("x*={%g, %g} rc=%d msg='%s' evals=%d\n", result.x[0], result.x[1], rc, result.message, result.nf);

c/examples/cobyla/cobyla_example.c

@@ -16,6 +16,17 @@ static void fun(const double x[], double *f, double constr[], const void *data)
   (void)data;
 }
 
+
+static void callback(const int n, const double x[], const double f, const int nf, const int tr, const double cstrv, const int m_nlcon, const double nlconstr[], bool *terminate)
+{
+  (void)n;
+  printf("progress: x=[%g;%g] f=%g cstrv=%g nf=%d tr=%d\n", x[0], x[1], f, cstrv, nf, tr);
+  *terminate = 0;
+  (void)m_nlcon;
+  (void)nlconstr;
+}
+
+
 int main(int argc, char * argv[])
 {
   (void)argc;
@@ -31,6 +42,7 @@ int main(int argc, char * argv[])
   options.iprint = PRIMA_MSG_EXIT;
   options.rhoend= 1e-3;
   options.maxfun = 200*n;
+  options.callback = &callback;
   problem.m_nlcon = M_NLCON;
   // x1<=4, x2<=3, x1+x2<=10
   problem.m_ineq = 3;

c/examples/lincoa/lincoa_example.c

@@ -12,6 +12,15 @@ static void fun(const double x[], double *f, const void *data)
   (void)data;
 }
 
+static void callback(int n, const double x[], double f, int nf, int tr, double cstrv, int m_nlcon, const double nlconstr[], bool *terminate)
+{
+  (void)n;
+  printf("progress: x=[%g;%g] f=%g cstrv=%g nf=%d tr=%d\n", x[0], x[1], f, cstrv, nf, tr);
+  (void)m_nlcon;
+  (void)nlconstr;
+  *terminate = 0;
+}
+
 int main(int argc, char * argv[])
 {
   (void)argc;
@@ -27,6 +36,7 @@ int main(int argc, char * argv[])
   options.iprint = PRIMA_MSG_EXIT;
   options.rhoend= 1e-3;
   options.maxfun = 200*n;
+  options.callback = &callback;
   // x1<=4, x2<=3, x1+x2<=10
   problem.m_ineq = 3;
   double Aineq[3*2] = {1.0, 0.0,

c/examples/newuoa/newuoa_example.c

@@ -12,6 +12,15 @@ static void fun(const double x[], double *f, const void *data)
   (void)data;
 }
 
+static void callback(int n, const double x[], double f, int nf, int tr, double cstrv, int m_nlcon, const double nlconstr[], bool *terminate)
+{
+  (void)n;
+  printf("progress: x=[%g;%g] f=%g cstrv=%g nf=%d tr=%d\n", x[0], x[1], f, cstrv, nf, tr);
+  (void)m_nlcon;
+  (void)nlconstr;
+  *terminate = 0;
+}
+
 int main(int argc, char * argv[])
 {
   (void)argc;
@@ -27,6 +36,7 @@ int main(int argc, char * argv[])
   options.iprint = PRIMA_MSG_EXIT;
   options.rhoend= 1e-3;
   options.maxfun = 200*n;
+  options.callback = &callback;
   prima_result_t result;
   const int rc = prima_minimize(PRIMA_NEWUOA, &problem, &options, &result);
   printf("x*={%g, %g} rc=%d msg='%s' evals=%d\n", result.x[0], result.x[1], rc, result.message, result.nf);

c/examples/uobyqa/uobyqa_example.c

@@ -12,6 +12,15 @@ static void fun(const double x[], double *f, const void *data)
   (void)data;
 }
 
+static void callback(int n, const double x[], double f, int nf, int tr, double cstrv, const int m_nlcon, const double nlconstr[], bool *terminate)
+{
+  (void)n;
+  printf("progress: x=[%g;%g] f=%g cstrv=%g nf=%d tr=%d\n", x[0], x[1], f, cstrv, nf, tr);
+  *terminate = 0;
+  (void)m_nlcon;
+  (void)nlconstr;
+}
+
 int main(int argc, char * argv[])
 {
   (void)argc;
@@ -27,6 +36,7 @@ int main(int argc, char * argv[])
   options.iprint = PRIMA_MSG_EXIT;
   options.rhoend= 1e-3;
   options.maxfun = 200*n;
+  options.callback = &callback;
   prima_result_t result;
   const int rc = prima_minimize(PRIMA_UOBYQA, &problem, &options, &result);
   printf("x*={%g, %g} rc=%d msg='%s' evals=%d\n", result.x[0], result.x[1], rc, result.message, result.nf);

c/include/prima/prima.h

@@ -3,6 +3,8 @@
 #ifndef PRIMA_H
 #define PRIMA_H
 
+#include <stdbool.h>
+
 #ifdef __cplusplus
 extern "C" {
 #endif
@@ -49,6 +51,7 @@ typedef enum
   PRIMA_NO_SPACE_BETWEEN_BOUNDS = 6,
   PRIMA_DAMAGING_ROUNDING = 7,
   PRIMA_ZERO_LINEAR_CONSTRAINT = 8,
+  PRIMA_USER_STOP = 9,
   PRIMA_INVALID_INPUT = 100,
   PRIMA_ASSERTION_FAILS = 101,
   PRIMA_VALIDATION_FAILS = 102,
@@ -67,7 +70,7 @@ PRIMAC_API
 const char *prima_get_rc_string(const prima_rc_t rc);
 
 /*
- * A function as required by solvers
+ * The objective function as required by solvers
  *
  * x     : on input, then vector of variables (should not be modified)
  * f     : on output, the value of the function
@@ -80,6 +83,21 @@ const char *prima_get_rc_string(const prima_rc_t rc);
 typedef void (*prima_obj_t)(const double x[], double *f, const void *data);
 typedef void (*prima_objcon_t)(const double x[], double *f, double constr[], const void *data);
 
+/* An optional callback function to report algorithm progress
+ *
+ * n     : number of variables
+ * x     : the current best point
+ * f     : the function value of the best point
+ * nf    : number of objective function calls
+ * tr    : iteration number
+ * cstrv : the constraint value verified by the current best point
+ * m_nlcon : number of non-linear constraints (cobyla only)
+ * nlconstr : non-linear constraints values verified by the current best point (cobyla only)
+ * terminate : a boolean to ask from early optimization exit
+*/
+typedef void (*prima_callback_t)(const int n, const double x[], const double f, int nf, int tr,
+                               const double cstrv, int m_nlcon, const double nlconstr[], bool *terminate);
+
 
 typedef struct {
 
@@ -105,6 +123,9 @@ typedef struct {
   // user-data, will be passed through the objective function callback
   void *data;
 
+  // callback function to report algorithm progress (default=NULL)
+  prima_callback_t callback;
+
 } prima_options_t;
 
 /* Initialize problem */