From c85e631a435a2b66a2bf30805f1bfd565da62eab Mon Sep 17 00:00:00 2001 From: zaikunzhang Date: Mon, 29 Jan 2024 13:09:36 +0800 Subject: [PATCH] 240129.130936.HKT improve stress. --- c/tests/data.c | 29 +++++----- c/tests/stress.c | 137 ++++++++++++++++++++++++++++------------------- 2 files changed, 100 insertions(+), 66 deletions(-) diff --git a/c/tests/data.c b/c/tests/data.c index f3c2972773..8754b8b2b6 100644 --- a/c/tests/data.c +++ b/c/tests/data.c @@ -5,6 +5,7 @@ #include #include +// Make PRIMA available #include "prima/prima.h" #define M_NLCON 1 @@ -14,6 +15,7 @@ int debug = 0; static int int_data = 0xff; void * data_ref = &int_data; +// Objective function for unconstrained, bound constrained, and linearly-constrained problems static void fun(const double x[], double *f, const void *data) { const double x1 = x[0]; @@ -35,6 +37,7 @@ static void fun(const double x[], double *f, const void *data) } } +// Objective & constraint function for nonlinearly-constrained problems static void fun_con(const double x[], double *f, double constr[], const void *data) { const double x1 = x[0]; @@ -57,6 +60,7 @@ static void fun_con(const double x[], double *f, double constr[], const void *da } } +// Main function int main(int argc, char * argv[]) { char *algo = "uobyqa"; @@ -67,8 +71,16 @@ int main(int argc, char * argv[]) if (argc > 2) debug = (strcmp(argv[2], "debug") == 0); - printf("debug = %d\n", debug); + printf("Debug = %d\n", debug); + // Define the options for the algorithm + prima_options_t options; + prima_init_options(&options); + options.iprint = PRIMA_MSG_RHO; + options.maxfun = 500*n; + options.data = data_ref; + + // Data for the problem double x0[] = {0.0, 0.0}; double xl[] = {-6.0, @@ -82,19 +94,10 @@ int main(int argc, char * argv[]) 3.0, 10.0}; + // Define the algorithm and the problem according to `algo` prima_problem_t problem; prima_init_problem(&problem, n); problem.x0 = x0; - - prima_options_t options; - prima_init_options(&options); - options.iprint = PRIMA_MSG_RHO; - options.maxfun = 500*n; - options.data = data_ref; - - prima_result_t result; - - // Define the algorithm and the problem according to `algo` if(strcmp(algo, "uobyqa") == 0) { algorithm = PRIMA_UOBYQA; @@ -135,11 +138,13 @@ int main(int argc, char * argv[]) } else { - printf("Invalid algorithm!\n"); + printf("Invalid algorithm %s!\n", algo); return 1; } // Call the solver + prima_result_t result; + int rc = prima_minimize(algorithm, &problem, &options, &result); // Print the result diff --git a/c/tests/stress.c b/c/tests/stress.c index fae6ce15b7..42a1a16eef 100644 --- a/c/tests/stress.c +++ b/c/tests/stress.c @@ -1,12 +1,14 @@ -// A stress test on excessively large problems. +// A stress test on excessively large problems -#include "prima/prima.h" -#include #include +#include #include #include #include +// Make PRIMA available +#include "prima/prima.h" + #define MIN(x, y) (((x) < (y)) ? (x) : (y)) #define N_MAX 2000 @@ -23,53 +25,58 @@ static double random_gen(double a, double b) return a + rand() * (b - a) / RAND_MAX; } +// Objective function for unconstrained, bound constrained, and linearly-constrained problems static void fun(const double x[], double *f, const void *data) { - // Rosenbrock function + // Objective: Rosenbrock function *f = 0.0; for (int i = 0; i < n-1; ++ i) *f += (x[i] - 1.0) * (x[i] - 1.0) + alpha * (x[i+1] - x[i]*x[i]) * (x[i+1] - x[i]*x[i]); - static int count = 0; + static int nf = 0; if (debug) { - ++ count; - printf("count=%d\n", count); + ++ nf; + printf("Number of function evaluations = %d\n", nf); } (void)data; } +// Objective & constraint function for nonlinearly-constrained problems static void fun_con(const double x[], double *f, double constr[], const void *data) { - // Rosenbrock function + // Objective: Rosenbrock function *f = 0.0; for (int i = 0; i < n-1; ++ i) *f += (x[i] - 1.0) * (x[i] - 1.0) + alpha * (x[i+1] - x[i]*x[i]) * (x[i+1] - x[i]*x[i]); - // x_{i+1} <= x_i^2 + + // Constraint: x_{i+1} <= x_i^2 for (int i = 0; i < MIN(M_NLCON, n-1); ++ i) constr[i] = x[i+1] - x[i] * x[i]; - static int count = 0; + static int nf = 0; if (debug) { - ++ count; - printf("count=%d\n", count); + ++ nf; + printf("Number of function evaluations = %d\n", nf); } (void)data; } +// Main function int main(int argc, char * argv[]) { - char *algo = "bobyqa"; + char *algo = "uobyqa"; + prima_algorithm_t algorithm = PRIMA_UOBYQA; if (argc > 1) algo = argv[1]; - printf("algo=%s\n", algo); + printf("Algorithm = %s\n", algo); if (argc > 2) debug = (strcmp(argv[2], "debug") == 0); - printf("debug=%d\n", debug); + printf("Debug = %d\n", debug); - // set seed to year/week + // Set the random seed to year/week char buf[10] = {0}; time_t t = time(NULL); struct tm *tmp = localtime(&t); @@ -77,75 +84,97 @@ int main(int argc, char * argv[]) if (!rc) return 1; unsigned seed = atoi(buf); - printf("seed=%d\n", seed); + printf("Random seed = %d\n", seed); srand(seed); - double x0[N_MAX]; - double xl[N_MAX]; - double xu[N_MAX]; - prima_problem_t problem; - prima_init_problem(&problem, N_MAX); - problem.x0 = x0; - problem.calcfc = &fun_con; - problem.calfun = &fun; + // Define the options for the algorithm prima_options_t options; prima_init_options(&options); options.iprint = PRIMA_MSG_RHO; options.maxfun = 500*N_MAX; + + // Data for the problem + double x0[N_MAX]; + double xl[N_MAX]; + double xu[N_MAX]; + for (int i = 0; i < N_MAX; ++ i) + { + x0[i] = random_gen(-1.0, 1.0); + xl[i] = random_gen(-2.0, -1.0); + xu[i] = random_gen(1.0, 2.0); + } double *Aineq = malloc(N_MAX*M_INEQ_MAX*sizeof(double)); + for (int i = 0; i < N_MAX; ++ i) + { + for (int j = 0; j < m_ineq; ++ j) + Aineq[j*N_MAX+i] = random_gen(-1.0, 1.0); + } double bineq[M_INEQ_MAX]; - problem.Aineq = Aineq; - problem.bineq = bineq; - problem.xl = xl; - problem.xu = xu; for (int j = 0; j < M_INEQ_MAX; ++ j) bineq[j] = random_gen(-1.0, 1.0); - for (int i = 0; i < N_MAX; ++ i) + // Define the algorithm and the problem according to `algo` + prima_problem_t problem; + prima_init_problem(&problem, N_MAX); + problem.x0 = x0; + if(strcmp(algo, "uobyqa") == 0) { - for (int j = 0; j < m_ineq; ++ j) - Aineq[j*N_MAX+i] = random_gen(-1.0, 1.0); - x0[i] = random_gen(-1.0, 1.0); - xl[i] = -1.0; - xu[i] = 1.0; + algorithm = PRIMA_UOBYQA; + problem.n = 100; + problem.calfun = &fun; } - prima_algorithm_t algorithm = 0; - prima_result_t result; - if(strcmp(algo, "bobyqa") == 0) + else if(strcmp(algo, "newuoa") == 0) { - algorithm = PRIMA_BOBYQA; + algorithm = PRIMA_NEWUOA; problem.n = 1600; + problem.calfun = &fun; } - else if(strcmp(algo, "cobyla") == 0) + else if(strcmp(algo, "bobyqa") == 0) { - algorithm = PRIMA_COBYLA; - problem.n = 800; - problem.m_nlcon = M_NLCON; - problem.m_ineq = 600; + algorithm = PRIMA_BOBYQA; + problem.n = 1600; + problem.calfun = &fun; + problem.xl = xl; + problem.xu = xu; } else if(strcmp(algo, "lincoa") == 0) { algorithm = PRIMA_LINCOA; problem.n = 1000; problem.m_ineq = 1000; + problem.calfun = &fun; + problem.xl = xl; + problem.xu = xu; + problem.Aineq = Aineq; + problem.bineq = bineq; } - else if(strcmp(algo, "newuoa") == 0) - { - algorithm = PRIMA_NEWUOA; - problem.n = 1600; - } - else if(strcmp(algo, "uobyqa") == 0) + else if(strcmp(algo, "cobyla") == 0) { - algorithm = PRIMA_UOBYQA; - problem.n = 100; + algorithm = PRIMA_COBYLA; + problem.n = 800; + problem.m_nlcon = M_NLCON; + problem.m_ineq = 600; + problem.calcfc = &fun_con; + problem.xl = xl; + problem.xu = xu; + problem.Aineq = Aineq; + problem.bineq = bineq; } else { - printf("incorrect algo\n"); + printf("Invalid algorithm %s!\n", algo); return 1; } + + // Call the solver + prima_result_t result; rc = prima_minimize(algorithm, &problem, &options, &result); - 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, result.message, result.nf); + + // Print the result + 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, result.message, result.nf); + + // Free the result prima_free_result(&result); + return 0; }