/* nag_opt_handle_disable (e04tcc) Example Program.
*
* Copyright 2021 Numerical Algorithms Group.
*
* Mark 27.3, 2021.
*/
#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
struct usr_data {
Integer nres;
double *t, *y;
};
int main(void) {
/* This example shows how to use the NAG Optimization Modeling Suite
* to edit a nonlinear least squares problem.
* In particular, the outlier residuals are disabled without having to
* redefine the full problem and comparing the solutions is easy.
*
* We try to fit the following model with 5 parameters:
* f(t) = at^2 + bt + c + d sin(omega t)
* to some noisy data (30 points) with 2 outliers (point 10 and 20)
* The data was generated using the values:
* a = 0.3
* b = 1.0
* c = 0.01
* d = 0.2
* omega = 5.0 */
Integer nvar, nres, isparse, nnzrd, idx[2];
double x[5];
double rinfo[100], stats[100];
double *rx;
void *handle = 0;
Integer exit_status = 0, i;
struct usr_data ud;
/* Nag Types */
Nag_Comm comm;
NagError fail;
printf("nag_opt_handle_disable (e04tcc) Example Program Results\n\n");
fflush(stdout);
exit_status = 0;
scanf(" %*[^\n]"); /* Skip heading in data file */
/* Read the data dimension */
scanf("%" NAG_IFMT " %*[^\n]", &nres);
/* Allocate the user data memory */
ud.nres = nres;
if (!(ud.t = NAG_ALLOC(nres, double)) || !(ud.y = NAG_ALLOC(nres, double)) ||
!(rx = NAG_ALLOC(nres, double))) {
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
/* Read the data points */
for (i = 0; i < nres; i++) {
scanf("%lf", &ud.t[i]);
}
scanf("%*[^\n]");
for (i = 0; i < nres; i++) {
scanf("%lf", &ud.y[i]);
}
/* nag_opt_handle_init (e04rac).
* Initialize the handle
*/
nvar = 5;
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 printing options
*/
nag_opt_handle_opt_set(handle, "Print Level = 1", NAGERR_DEFAULT);
nag_opt_handle_opt_set(handle, "Print Options = No", NAGERR_DEFAULT);
nag_opt_handle_opt_set(handle, "Print Solution = X", NAGERR_DEFAULT);
printf(" First solve the problem with the outliers\n\n");
printf(" --------------------------------------------------------\n");
fflush(stdout);
/* nag_opt_handle_solve_bxnl (e04ggc)
* Call the solver
*/
comm.p = &ud;
for (i = 0; i < nvar; i++) {
x[i] = 1.0;
}
INIT_FAIL(fail);
nag_opt_handle_solve_bxnl(handle, lsqfun, lsqgrd, NULLFN, NULLFN, NULLFN,
nvar, x, nres, rx, rinfo, stats, &comm, &fail);
printf(" --------------------------------------------------------\n");
printf("\n Now remove the outlier residuals from the problem handle\n\n");
printf(" --------------------------------------------------------\n");
fflush(stdout);
/* nag_opt_handle_disable (e04tcc)
* Disable the two outlier residuals */
idx[0] = 10;
idx[1] = 20;
nag_opt_handle_disable(handle, "NLS", 2, idx, NAGERR_DEFAULT);
/* Solve the problem again */
for (i = 0; i < nvar; i++) {
x[i] = 1.0;
}
INIT_FAIL(fail);
nag_opt_handle_solve_bxnl(handle, lsqfun, lsqgrd, NULLFN, NULLFN, NULLFN,
nvar, x, nres, rx, rinfo, stats, &comm, &fail);
printf(" --------------------------------------------------------\n");
printf("\n Assuming the outliers points are measured again\n");
printf(" we can enable the residuals and adjust the values\n\n");
printf(" --------------------------------------------------------\n");
fflush(stdout);
/* Fix the first variable to its known value of 0.3
* enable the residuals and adjust the values in the data */
nag_opt_handle_set_bound(handle, "variable", 1, 0.3, 0.3, NAGERR_DEFAULT);
nag_opt_handle_enable(handle, "NLS", 2, idx, NAGERR_DEFAULT);
ud.y[9] = -0.515629;
ud.y[19] = 0.54920;
/* Solve the problem again */
for (i = 0; i < nvar; i++) {
x[i] = 1.0;
}
INIT_FAIL(fail);
nag_opt_handle_solve_bxnl(handle, lsqfun, lsqgrd, NULLFN, NULLFN, NULLFN,
nvar, x, nres, rx, rinfo, stats, &comm, &fail);
END:
NAG_FREE(ud.t);
NAG_FREE(ud.y);
NAG_FREE(rx);
/* nag_opt_handle_free (e04rzc).
* Destroy the problem handle and deallocate all the memory. */
if (handle)
nag_opt_handle_free(&handle, NAGERR_DEFAULT);
return exit_status;
}
static void NAG_CALL lsqfun(Integer nvar, const double x[], Integer nres,
double rx[], Integer *inform, Nag_Comm *comm) {
Integer i;
struct usr_data *ud;
/* Interrupt the solver if the comm structure is not correctly initialized */
if (!comm || !(comm->p)) {
*inform = -1;
return;
}
*inform = 0;
ud = (struct usr_data *)comm->p;
/* Fill the residuals values */
for (i = 0; i < nres; i++) {
rx[i] = x[0] * pow(ud->t[i], 2) + x[1] * ud->t[i] + x[2] +
x[3] * sin(x[4] * ud->t[i]) - ud->y[i];
}
}
static void NAG_CALL lsqgrd(Integer nvar, const double x[], Integer nres,
Integer nnzrd, double rdx[], Integer *inform,
Nag_Comm *comm) {
Integer i;
struct usr_data *ud;
/* Interrupt the solver if the comm structure is not correctly initialized */
if (!comm || !(comm->p)) {
*inform = -1;
return;
}
*inform = 0;
ud = (struct usr_data *)comm->p;
for (i = 0; i < nres; i++) {
rdx[i * nvar + 0] = pow(ud->t[i], 2);
rdx[i * nvar + 1] = ud->t[i];
rdx[i * nvar + 2] = 1.0;
rdx[i * nvar + 3] = sin(x[4] * ud->t[i]);
rdx[i * nvar + 4] = x[3] * ud->t[i] * cos(x[4] * ud->t[i]);
}
}