/* nag_opt_handle_solve_nldf (e04gnc) Example Program.
*
* Copyright 2023 Numerical Algorithms Group.
*
* Mark 29.1, 2023.
*/
#include <math.h>
#include <nag.h>
#include <stdio.h>
#ifdef __cplusplus
extern "C" {
#endif
static void NAG_CALL lsqfun(Integer nvar, const double x[], Integer nres,
double rx[], Integer *inform, Nag_Comm *comm);
static void NAG_CALL lsqgrd(Integer nvar, const double x[], Integer nres,
Integer nnzrd, double rdx[], Integer *inform,
Nag_Comm *comm);
#ifdef __cplusplus
}
#endif
/* Model data
* y = f(t;A,w) := A*t*sin(-w*t)
* Data generated using A = 0.1 and w = 0.8 and added random noise
* Residual 4, 12, 16 and 20 are outliers
*/
const double t[24] = {1.0E+00, 2.0E+00, 3.0E+00, 4.0E+00, 5.0E+00, 6.0E+00,
7.0E+00, 8.0E+00, 9.0E+00, 1.0E+01, 1.1E+01, 1.2E+01,
1.3E+01, 1.4E+01, 1.5E+01, 1.6E+01, 1.7E+01, 1.8E+01,
1.9E+01, 2.0E+01, 2.1E+01, 2.2E+01, 2.3E+01, 2.4E+01};
const double y[24] = { 0.0523, -0.1442, -0.0422, 1.8106, 0.3271, 0.4684,
0.4593, -0.0169, -0.7811, -1.1356, -0.5343, -3.0043,
1.1832, 1.5153, 0.7120, -2.2923, -1.4871, -1.7083,
-0.9936, -5.2873, 1.7555, 2.0642, 0.9499, -0.6234};
int main(void) {
const double infbnd = 1.0e20;
Integer nvar, nres, nclin, nnza, isparse, nnzrd, x_idx, idlc;
Integer irowa[2] = {1, 1};
Integer icola[2] = {1, 2};
double x[2] = {0.3, 0.7};
double rinfo[100], stats[100];
double blx[2] = {-1.0, 0.0};
double bux[2] = {infbnd, 1.0};
double bla[1] = {0.2};
double bua[1] = {1.0};
double a[2] = {1.0, 1.0};
double *rx = 0;
void *handle = 0;
Integer exit_status = 0;
/* Nag Types */
Nag_Comm comm;
NagError fail;
printf("nag_opt_handle_solve_nldf (e04gnc) Example Program Results\n\n");
fflush(stdout);
nvar = 2;
nres = 24;
/* Allocate memory */
if (!(comm.user = NAG_ALLOC(2 * nres, double)) ||
!(rx = NAG_ALLOC(nres, double))) {
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
/* Fill the problem data structure
* copy t[]
*/
memcpy(&comm.user[0], t, 24 * sizeof(double));
/* copy y[] = */
memcpy(&comm.user[24], y, 24 * sizeof(double));
/* nag_opt_handle_init (e04rac).
* Initialize the handle
*/
nag_opt_handle_init(&handle, nvar, NAGERR_DEFAULT);
/* nag_opt_handle_set_nlnls (e04rmc)
* Define residuals structure, isparse=0 means the residual structure is
* dense => irowrd and icolrd arguments can be NULL
*/
isparse = 0;
nnzrd = 1;
nag_opt_handle_set_nlnls(handle, nres, isparse, nnzrd, NULL, NULL,
NAGERR_DEFAULT);
/* nag_opt_handle_opt_set (e04zmc)
* Set options
*/
/* Use least-squares loss function */
nag_opt_handle_opt_set(handle, "NLDF Loss Function Type = L2",
NAGERR_DEFAULT);
/* Reduce printed output info */
nag_opt_handle_opt_set(handle, "Print Level = 1", NAGERR_DEFAULT);
nag_opt_handle_opt_set(handle, "Print Options = No", NAGERR_DEFAULT);
/* Define bounds */
/* nag_opt_handle_set_simplebounds (e04rhc) */
nag_opt_handle_set_simplebounds(handle, nvar, blx, bux, NAGERR_DEFAULT);
printf(" First solve the problem using least squares loss function\n");
printf(" --------------------------------------------------------\n");
fflush(stdout);
/* nag_opt_handle_set_linconstr (e04rjc)
* Define linear constraints */
nclin = 1;
nnza = 2;
idlc = 0;
nag_opt_handle_set_linconstr(handle, nclin, bla, bua, nnza, irowa, icola, a,
&idlc, NAGERR_DEFAULT);
/* nag_opt_handle_solve_nldf (e04gnc)
* Call the solver
*/
INIT_FAIL(fail);
nag_opt_handle_solve_nldf(handle, lsqfun, lsqgrd, NULLFN, NULLFN, NULLFN,
nvar, x, nres, rx, rinfo, stats, &comm, &fail);
if (fail.code != NE_NOERROR && fail.code != NW_NOT_CONVERGED) {
printf("Error from nag_opt_handle_solve_nldf (e04gnc).\n%s\n",
fail.message);
exit_status = 1;
} else {
/* Print the solution point */
printf("\n Optimal parameters X (Least Squares):\n");
printf(" x_idx Value\n");
for (x_idx = 0; x_idx < nvar; x_idx++) {
printf(" %5" NAG_IFMT " %7.2E\n", x_idx + 1, x[x_idx]);
}
}
printf("\n Now solve the problem using SmoothL1 loss function\n");
printf(" --------------------------------------------------------\n");
fflush(stdout);
/* nag_opt_handle_opt_set (e04zmc)
* Set options
*/
/* Use SmoothL1 loss function */
nag_opt_handle_opt_set(handle, "NLDF Loss Function Type = SmoothL1",
NAGERR_DEFAULT);
/* nag_opt_handle_solve_nldf (e04gnc)
* Call the solver
*/
INIT_FAIL(fail);
x[0] = 0.3;
x[1] = 0.7;
nag_opt_handle_solve_nldf(handle, lsqfun, lsqgrd, NULLFN, NULLFN, NULLFN,
nvar, x, nres, rx, rinfo, stats, &comm, &fail);
if (fail.code != NE_NOERROR && fail.code != NW_NOT_CONVERGED) {
printf("Error from nag_opt_handle_solve_nldf (e04gnc).\n%s\n",
fail.message);
exit_status = 1;
} else {
/* Print the solution point */
/* Print the solution point */
printf("\n Optimal parameters X (SmoothL1):\n");
printf(" x_idx Value\n");
for (x_idx = 0; x_idx < nvar; x_idx++) {
printf(" %5" NAG_IFMT " %7.2E\n", x_idx + 1, x[x_idx]);
}
}
END:
/* Free allocated memory */
if (handle)
/* nag_opt_handle_free (e04rzc).
* Free the problem handle and deallocate all the memory used
*/
nag_opt_handle_free(&handle, NAGERR_DEFAULT);
if (comm.user)
NAG_FREE(comm.user);
NAG_FREE(rx);
return exit_status;
}
static void NAG_CALL lsqfun(Integer nvar, const double x[], Integer nres,
double rx[], Integer *inform, Nag_Comm *comm) {
int i;
/* Interrupt the solver if the comm structure is not correctly initialized */
if (!comm || !(comm->user)) {
*inform = -1;
return;
}
/* Interrupt the solver if the data does not correspond to the problem */
if (nvar != 2 || nres != 24) {
*inform = -1;
return;
}
*inform = 0;
/* Fill the residuals values */
for (i = 0; i < nres; i++) {
rx[i] = comm->user[nres + i];
rx[i] = rx[i] - x[0] * comm->user[i] * sin(-x[1] * comm->user[i]);
}
}
static void NAG_CALL lsqgrd(Integer nvar, const double x[], Integer nres,
Integer nnzrd, double rdx[], Integer *inform,
Nag_Comm *comm) {
int i;
/* Interrupt the solver if the comm structure is not correctly initialized */
if (!comm || !(comm->user)) {
*inform = -1;
return;
}
/* Interrupt the solver if the data does not correspond to the problem */
if (nvar != 2 || nres != 24 || nnzrd != nres * nvar) {
*inform = -1;
return;
}
*inform = 0;
for (i = 0; i < nres; i++) {
rdx[i * nvar + 0] = -comm->user[i] * sin(-x[1] * comm->user[i]);
rdx[i * nvar + 1] = comm->user[i] * comm->user[i] * x[0]
* cos(x[1] * comm->user[i]);
}
}