From 870793ed83413bd83e178fec148971fc18444272 Mon Sep 17 00:00:00 2001
From: Julien Schueller <schueller@phimeca.com>
Date: Tue, 24 Oct 2023 10:02:25 +0200
Subject: [PATCH] Add prima_result

---
 c/examples/bobyqa/bobyqa_example.c |  8 ++--
 c/examples/cobyla/cobyla_example.c | 10 ++--
 c/examples/lincoa/lincoa_example.c |  6 +--
 c/examples/newuoa/newuoa_example.c |  7 ++-
 c/examples/uobyqa/uobyqa_example.c |  7 ++-
 c/include/prima/prima.h            | 49 ++++++++++++-------
 c/prima.c                          | 76 +++++++++++++++++++++---------
 c/tests/data.c                     | 20 ++++----
 c/tests/stress.c                   | 23 ++++-----
 9 files changed, 117 insertions(+), 89 deletions(-)

diff --git a/c/examples/bobyqa/bobyqa_example.c b/c/examples/bobyqa/bobyqa_example.c
index ac438b76a9..f1c3578daf 100644
--- a/c/examples/bobyqa/bobyqa_example.c
+++ b/c/examples/bobyqa/bobyqa_example.c
@@ -18,16 +18,14 @@ int main(int argc, char * argv[])
   (void)argv;
   const int n = 2;
   double x[2] = {0.0, 0.0};
-  double f = 0.0;
-
   prima_options options;
   prima_init_options(&options);
   options.iprint = PRIMA_MSG_EXIT;
   options.rhoend= 1e-3;
   options.maxfun = 200*n;
-  int nf = 0;
-  const int rc = prima_bobyqa(&fun, n, x, &f, &nf, &options);
+  prima_result result;
+  const int rc = prima_bobyqa(&fun, n, x, &options, &result);
   const char *msg = prima_get_rc_string(rc);
-  printf("x*={%g, %g} rc=%d msg='%s' evals=%d\n", x[0], x[1], rc, msg, nf);
+  printf("x*={%g, %g} rc=%d msg='%s' evals=%d\n", x[0], x[1], rc, msg, result.nf);
   return (fabs(x[0]-3)>2e-2 || fabs(x[1]-2)>2e-2);
 }
diff --git a/c/examples/cobyla/cobyla_example.c b/c/examples/cobyla/cobyla_example.c
index 21bc5fa280..340b5c0c3b 100644
--- a/c/examples/cobyla/cobyla_example.c
+++ b/c/examples/cobyla/cobyla_example.c
@@ -21,9 +21,6 @@ int main(int argc, char * argv[])
   const int m_nlcon = 1;
   const int n = 2;
   double x[2] = {0.0, 0.0};
-  double f = 0.0;
-  double cstrv = 0.0;
-  double nlconstr[m_nlcon];
   prima_options options;
   prima_init_options(&options);
   options.iprint = PRIMA_MSG_EXIT;
@@ -44,9 +41,10 @@ int main(int argc, char * argv[])
   double xu[2] = {6.0, 6.0};
   options.xl = xl;
   options.xu = xu;
-  int nf = 0;
-  const int rc = prima_cobyla(&fun, n, x, &f, &cstrv, nlconstr, &nf, &options);
+  prima_result result;
+  const int rc = prima_cobyla(&fun, n, x, &options, &result);
   const char *msg = prima_get_rc_string(rc);
-  printf("x*={%g, %g} f*=%g cstrv=%g nlconstr=%g rc=%d msg='%s' evals=%d\n", x[0], x[1], f, cstrv, nlconstr[0], rc, msg, nf);
+  printf("x*={%g, %g} f*=%g cstrv=%g nlconstr=%g rc=%d msg='%s' evals=%d\n", x[0], x[1], result.f, result.cstrv, result.nlconstr[0], rc, msg, result.nf);
+  prima_free_result(&result);
   return (fabs(x[0]-3)>2e-2 || fabs(x[1]-2)>2e-2);
 }
diff --git a/c/examples/lincoa/lincoa_example.c b/c/examples/lincoa/lincoa_example.c
index e0f0ae9586..f800f8f4f0 100644
--- a/c/examples/lincoa/lincoa_example.c
+++ b/c/examples/lincoa/lincoa_example.c
@@ -39,9 +39,9 @@ int main(int argc, char * argv[])
   double xu[2] = {6.0, 6.0};
   options.xl = xl;
   options.xu = xu;
-  int nf = 0;
-  const int rc = prima_lincoa(&fun, n, x, &f, &cstrv, &nf, &options);
+  prima_result result;
+  const int rc = prima_lincoa(&fun, n, x, &options, &result);
   const char *msg = prima_get_rc_string(rc);
-  printf("x*={%g, %g} f*=%g cstrv=%g rc=%d msg='%s' evals=%d\n", x[0], x[1], f, cstrv, rc, msg, nf);
+  printf("x*={%g, %g} f*=%g cstrv=%g rc=%d msg='%s' evals=%d\n", x[0], x[1], result.f, result.cstrv, rc, msg, result.nf);
   return (fabs(x[0]-3)>2e-2 || fabs(x[1]-2)>2e-2);
 }
diff --git a/c/examples/newuoa/newuoa_example.c b/c/examples/newuoa/newuoa_example.c
index 4d65f75880..9e1d7d72ca 100644
--- a/c/examples/newuoa/newuoa_example.c
+++ b/c/examples/newuoa/newuoa_example.c
@@ -18,15 +18,14 @@ int main(int argc, char * argv[])
   (void)argv;
   const int n = 2;
   double x[2] = {0.0, 0.0};
-  double f = 0.0;
   prima_options options;
   prima_init_options(&options);
   options.iprint = PRIMA_MSG_EXIT;
   options.rhoend= 1e-3;
   options.maxfun = 200*n;
-  int nf = 0;
-  const int rc = prima_newuoa(&fun, n, x, &f, &nf, &options);
+  prima_result result;
+  const int rc = prima_newuoa(&fun, n, x, &options, &result);
   const char *msg = prima_get_rc_string(rc);
-  printf("x*={%g, %g} rc=%d msg='%s' evals=%d\n", x[0], x[1], rc, msg, nf);
+  printf("x*={%g, %g} rc=%d msg='%s' evals=%d\n", x[0], x[1], rc, msg, result.nf);
   return (fabs(x[0]-3)>2e-2 || fabs(x[1]-2)>2e-2);
 }
diff --git a/c/examples/uobyqa/uobyqa_example.c b/c/examples/uobyqa/uobyqa_example.c
index 868bc12bcd..066f26992d 100644
--- a/c/examples/uobyqa/uobyqa_example.c
+++ b/c/examples/uobyqa/uobyqa_example.c
@@ -18,15 +18,14 @@ int main(int argc, char * argv[])
   (void)argv;
   const int n = 2;
   double x[2] = {0.0, 0.0};
-  double f = 0.0;
   prima_options options;
   prima_init_options(&options);
   options.iprint = PRIMA_MSG_EXIT;
   options.rhoend= 1e-3;
   options.maxfun = 200*n;
-  int nf = 0;
-  const int rc = prima_uobyqa(&fun, n, x, &f, &nf, &options);
+  prima_result result;
+  const int rc = prima_uobyqa(&fun, n, x, &options, &result);
   const char *msg = prima_get_rc_string(rc);
-  printf("x*={%g, %g} rc=%d msg='%s' evals=%d\n", x[0], x[1], rc, msg, nf);
+  printf("x*={%g, %g} rc=%d msg='%s' evals=%d\n", x[0], x[1], rc, msg, result.nf);
   return (fabs(x[0]-3)>2e-2 || fabs(x[1]-2)>2e-2);
 }
diff --git a/c/include/prima/prima.h b/c/include/prima/prima.h
index 0439844f48..22ea090fe2 100644
--- a/c/include/prima/prima.h
+++ b/c/include/prima/prima.h
@@ -135,46 +135,59 @@ typedef struct {
   
 } prima_options;
 
+
 /* Initialize option data */
 PRIMAC_API
 int prima_init_options(prima_options * options);
 
+
+typedef struct {
+  // objective value
+  double f;
+
+  // number of objective function calls
+  int nf;
+
+  // constraint violation (cobyla & lincoa)
+  double cstrv;
+
+  // non-linear constraint values, of size m_nlcon (cobyla only)
+  double *nlconstr;
+
+  // whether prima had to allocate nlconstr0 (private, do not use)
+  int _allocated_nlconstr;
+
+} prima_result;
+
+
+PRIMAC_API
+int prima_free_result(prima_result * result);
+
 /*
  * calfun    : function to minimize (see prima_obj)
+ * calcfc    : function to minimize and constraints (see prima_objcon)
  * n         : number of variables (>=0)
  * x         : on input, initial estimate
  *             on output, the solution
- * f         : objective value (output)
- * nf        : number of objective function calls (output)
  * options   : optimization options (see prima_options)
- * calcfc    : function to minimize and constraints (see prima_objcon)
- * cstrv     : constraint violation (output)
- * nlconstr  : non-linear constraint values of size m_nlcon (output)
- *
+ * result    : optimization results (see prima_result)
  * return    : see prima_rc enum for return codes
  */
 
 PRIMAC_API
-int prima_bobyqa(const prima_obj calfun, const int n, double x[], double *f,
-                 int *nf, prima_options * options);
+int prima_bobyqa(const prima_obj calfun, const int n, double x[], prima_options * options, prima_result * result);
 
 PRIMAC_API
-int prima_newuoa(const prima_obj calfun, const int n, double x[], double *f,
-                 int *nf, prima_options * options);
+int prima_newuoa(const prima_obj calfun, const int n, double x[], prima_options * options, prima_result * result);
 
 PRIMAC_API
-int prima_uobyqa(const prima_obj calfun, const int n, double x[], double *f,
-                 int *nf, prima_options * options);
+int prima_uobyqa(const prima_obj calfun, const int n, double x[], prima_options * options, prima_result * result);
 
 PRIMAC_API
-int prima_cobyla(const prima_objcon calcfc, const int n, double x[], double *f,
-                 double *cstrv, double nlconstr[],
-                 int *nf, prima_options * options);
+int prima_cobyla(const prima_objcon calcfc, const int n, double x[], prima_options * options, prima_result * result);
 
 PRIMAC_API
-int prima_lincoa(const prima_obj calfun, const int n, double x[], double *f,
-                 double *cstrv,
-                 int *nf, prima_options * options);
+int prima_lincoa(const prima_obj calfun, const int n, double x[], prima_options * options, prima_result * result);
 
 #ifdef __cplusplus
 }
diff --git a/c/prima.c b/c/prima.c
index d0813cabb7..87effc7faf 100644
--- a/c/prima.c
+++ b/c/prima.c
@@ -110,14 +110,43 @@ void prima_free_options(prima_options * opt)
   }
 }
 
+int prima_init_result(prima_result * result)
+{
+  result->f = 0.0;
+  result->nf = 0;
+  result->cstrv = 0.0;
+  result->nlconstr = NULL;
+  result->_allocated_nlconstr = 0;
+  return 0;
+}
+
+int prima_free_result(prima_result * result)
+{
+  if (result)
+  {
+    if (result->_allocated_nlconstr)
+    {
+      free(result->nlconstr);
+      result->nlconstr = NULL;
+      result->_allocated_nlconstr = 0;
+    }
+  }
+  return 0;
+}
+
 /* these functions just call the fortran compatibility layer and return the status code */
-int prima_cobyla(const prima_objcon calcfc, const int n, double x[], double *f,
-                 double *cstrv, double nlconstr[],
-                 int *nf, prima_options * opt)
+int prima_cobyla(const prima_objcon calcfc, const int n, double x[], prima_options * opt, prima_result * result)
 {
   int info = prima_check_options(opt, n, 1);
-  if (info == 0)
+  if ((info == 0) && result)
   {
+    prima_init_result(result);
+    // allocate nlconstr
+    result->nlconstr = (double*)calloc(opt->m_nlcon, sizeof(double));
+    if (!result->nlconstr)
+       return PRIMA_MEMORY_ALLOCATION_FAILS;
+    result->_allocated_nlconstr = 1;
+
     // evaluate f0, nlconstr0 if either one is not provided
     if (opt->f0 == NAN || !opt->nlconstr0)
     {
@@ -130,59 +159,60 @@ int prima_cobyla(const prima_objcon calcfc, const int n, double x[], double *f,
       }
       calcfc(x, &(opt->f0), opt->nlconstr0, opt->data);
     }
-    cobyla_c(opt->m_nlcon, calcfc, opt->data, n, x, f, cstrv, nlconstr,
+    cobyla_c(opt->m_nlcon, calcfc, opt->data, n, x, &(result->f), &(result->cstrv), result->nlconstr,
             opt->m_ineq, opt->Aineq, opt->bineq, opt->m_eq, opt->Aeq, opt->beq,
-            opt->xl, opt->xu, opt->f0, opt->nlconstr0, nf, opt->rhobeg, opt->rhoend, opt->ftarget, opt->maxfun, opt->iprint, &info);
+            opt->xl, opt->xu, opt->f0, opt->nlconstr0, &(result->nf), opt->rhobeg, opt->rhoend, opt->ftarget, opt->maxfun, opt->iprint, &info);
     prima_free_options(opt);
   }
   return info;
 }
 
-int prima_bobyqa(const prima_obj calfun, const int n, double x[], double *f,
-                 int *nf, prima_options * opt)
+int prima_bobyqa(const prima_obj calfun, const int n, double x[], prima_options * opt, prima_result * result)
 {
   int info = prima_check_options(opt, n, 1);
-  if (info == 0)
+  if ((info == 0) && result)
   {
-    bobyqa_c(calfun, opt->data, n, x, f, opt->xl, opt->xu, nf, opt->rhobeg, opt->rhoend, opt->ftarget, opt->maxfun, opt->npt < 0 ? 2*n+1 : opt->npt, opt->iprint, &info);
+    prima_init_result(result);
+    bobyqa_c(calfun, opt->data, n, x, &(result->f), opt->xl, opt->xu, &(result->nf), opt->rhobeg, opt->rhoend, opt->ftarget, opt->maxfun, opt->npt, opt->iprint, &info);
+    
     prima_free_options(opt);
   }
   return info;
 }
 
-int prima_newuoa(const prima_obj calfun, const int n, double x[], double *f,
-                 int *nf, prima_options * opt)
+int prima_newuoa(const prima_obj calfun, const int n, double x[], prima_options * opt, prima_result * result)
 {
   int info = prima_check_options(opt, n, 0);
-  if (info == 0)
+  if ((info == 0) && result)
   {
-    newuoa_c(calfun, opt->data, n, x, f, nf, opt->rhobeg, opt->rhoend, opt->ftarget, opt->maxfun, opt->npt < 0 ? 2*n+1 : opt->npt, opt->iprint, &info);
+    prima_init_result(result);
+    newuoa_c(calfun, opt->data, n, x, &(result->f), &(result->nf), opt->rhobeg, opt->rhoend, opt->ftarget, opt->maxfun, opt->npt, opt->iprint, &info);
     prima_free_options(opt);
   }
   return info;
 }
 
-int prima_uobyqa(const prima_obj calfun, const int n, double x[], double *f,
-                 int *nf, prima_options * opt)
+int prima_uobyqa(const prima_obj calfun, const int n, double x[], prima_options * opt, prima_result * result)
 {
   int info = prima_check_options(opt, n, 0);
-  if (info == 0)
+  if ((info == 0) && result)
   {
-    uobyqa_c(calfun, opt->data, n, x, f, nf, opt->rhobeg, opt->rhoend, opt->ftarget, opt->maxfun, opt->iprint, &info);
+    prima_init_result(result);
+    uobyqa_c(calfun, opt->data, n, x, &(result->f), &(result->nf), opt->rhobeg, opt->rhoend, opt->ftarget, opt->maxfun, opt->iprint, &info);
     prima_free_options(opt);
   }
   return info;
 }
 
-int prima_lincoa(const prima_obj calfun, const int n, double x[], double *f, double *cstrv,
-                 int *nf, prima_options * opt)
+int prima_lincoa(const prima_obj calfun, const int n, double x[], prima_options * opt, prima_result * result)
 {
   int info = prima_check_options(opt, n, 1);
-  if (info == 0)
+  if ((info == 0) && result)
   {
-    lincoa_c(calfun, opt->data, n, x, f, cstrv,
+    prima_init_result(result);
+    lincoa_c(calfun, opt->data, n, x, &(result->f), &(result->cstrv),
             opt->m_ineq, opt->Aineq, opt->bineq, opt->m_eq, opt->Aeq, opt->beq,
-            opt->xl, opt->xu, nf, opt->rhobeg, opt->rhoend, opt->ftarget, opt->maxfun, opt->npt < 0 ? 2*n+1 : opt->npt, opt->iprint, &info);
+            opt->xl, opt->xu, &(result->nf), opt->rhobeg, opt->rhoend, opt->ftarget, opt->maxfun, opt->npt, opt->iprint, &info);
     prima_free_options(opt);
   }
   return info;
diff --git a/c/tests/data.c b/c/tests/data.c
index 88a6840e8c..12e72ba3de 100644
--- a/c/tests/data.c
+++ b/c/tests/data.c
@@ -69,11 +69,6 @@ int main(int argc, char * argv[])
   double x[] = {0, 0};
   double xl[] = {-6.0, -6.0};
   double xu[] = {6.0, 6.0};
-  double f = 0.0;
-  double cstrv = 0.0;
-  double nlconstr[m_nlcon];
-  for (int j = 0; j < m_nlcon; ++ j)
-    nlconstr[j] = 0.0;
   prima_options options;
   prima_init_options(&options);
   options.iprint = PRIMA_MSG_RHO;
@@ -90,28 +85,28 @@ int main(int argc, char * argv[])
   options.bineq = bineq;
   options.xl = xl;
   options.xu = xu;
-  int nf = 0;
+  prima_result result;
   int rc = 0;
   if(strcmp(algo, "bobyqa") == 0)
   {
-    rc = prima_bobyqa(&fun, n, x, &f, &nf, &options);
+    rc = prima_bobyqa(&fun, n, x, &options, &result);
   }
   else if(strcmp(algo, "cobyla") == 0)
   {
     options.m_nlcon = m_nlcon;
-    rc = prima_cobyla(&fun_con, n, x, &f, &cstrv, nlconstr, &nf, &options);
+    rc = prima_cobyla(&fun_con, n, x, &options, &result);
   }
   else if(strcmp(algo, "lincoa") == 0)
   {
-    rc = prima_lincoa(&fun, n, x, &f, &cstrv, &nf, &options);
+    rc = prima_lincoa(&fun, n, x, &options, &result);
   }
   else if(strcmp(algo, "newuoa") == 0)
   {
-    rc = prima_newuoa(&fun, n, x, &f, &nf, &options);
+    rc = prima_newuoa(&fun, n, x, &options, &result);
   }
   else if(strcmp(algo, "uobyqa") == 0)
   {
-    rc = prima_uobyqa(&fun, n, x, &f, &nf, &options);
+    rc = prima_uobyqa(&fun, n, x, &options, &result);
   }
   else
   {
@@ -120,6 +115,7 @@ int main(int argc, char * argv[])
   }
   const char *msg = prima_get_rc_string(rc);
 
-  printf("f*=%g cstrv=%g nlconstr=%g rc=%d msg='%s' evals=%d\n", f, cstrv, nlconstr[0], rc, msg, nf);
+  printf("f*=%g cstrv=%g nlconstr=%g rc=%d msg='%s' evals=%d\n", result.f, result.cstrv, result.nlconstr ? result.nlconstr[0] : 0.0, rc, msg, result.nf);
+  prima_free_result(&result);
   return (fabs(x[0]-3)>2e-2 || fabs(x[1]-2)>2e-2);
 }
diff --git a/c/tests/stress.c b/c/tests/stress.c
index 9cc306dd08..55f180888d 100644
--- a/c/tests/stress.c
+++ b/c/tests/stress.c
@@ -82,10 +82,7 @@ int main(int argc, char * argv[])
   double x[n_max];
   double xl[n_max];
   double xu[n_max];
-  double f = 0.0;
-  double cstrv = 0.0;
-  double nlconstr[m_nlcon];
-    prima_options options;
+  prima_options options;
   prima_init_options(&options);
   options.iprint = PRIMA_MSG_RHO;
   options.maxfun = 500*n_max;
@@ -97,8 +94,6 @@ int main(int argc, char * argv[])
   options.xu = xu;
   for (int j = 0; j < m_ineq_max; ++ j)
     bineq[j] = random_gen(-1.0, 1.0);
-  for (int j = 0; j < m_nlcon; ++ j)
-    nlconstr[j] = 0.0;
 
   for (int i = 0; i < n_max; ++ i)
   {
@@ -108,35 +103,34 @@ int main(int argc, char * argv[])
     xl[i] = -1.0;
     xu[i] = 1.0;
   }
-  
-  int nf = 0;
+  prima_result result;
   if(strcmp(algo, "bobyqa") == 0)
   {
     n = 1600;
-    rc = prima_bobyqa(&fun, n, x, &f, &nf, &options);
+    rc = prima_bobyqa(&fun, n, x, &options, &result);
   }
   else if(strcmp(algo, "cobyla") == 0)
   {
     n = 800;
     options.m_nlcon = m_nlcon;
     options.m_ineq = 600;
-    rc = prima_cobyla(&fun_con, n, x, &f, &cstrv, nlconstr, &nf, &options);
+    rc = prima_cobyla(&fun_con, n, x, &options, &result);
   }
   else if(strcmp(algo, "lincoa") == 0)
   {
     n = 1000;
     options.m_ineq = 1000;
-    rc = prima_lincoa(&fun, n, x, &f, &cstrv, &nf, &options);
+    rc = prima_lincoa(&fun, n, x, &options, &result);
   }
   else if(strcmp(algo, "newuoa") == 0)
   {
     n = 1600;
-    rc = prima_newuoa(&fun, n, x, &f, &nf, &options);
+    rc = prima_newuoa(&fun, n, x, &options, &result);
   }
   else if(strcmp(algo, "uobyqa") == 0)
   {
     n = 100;
-    rc = prima_uobyqa(&fun, n, x, &f, &nf, &options);
+    rc = prima_uobyqa(&fun, n, x, &options, &result);
   }
   else
   {
@@ -145,6 +139,7 @@ int main(int argc, char * argv[])
   }
   const char *msg = prima_get_rc_string(rc);
 
-  printf("f*=%g cstrv=%g nlconstr=%g rc=%d msg='%s' evals=%d\n", f, cstrv, nlconstr[0], rc, msg, nf);
+  printf("f*=%g cstrv=%g nlconstr=%g rc=%d msg='%s' evals=%d\n", result.f, result.cstrv, result.nlconstr ? result.nlconstr[0] : 0.0, rc, msg, result.nf);
+  prima_free_result(&result);
   return 0;
 }