+/***********************************************************************/
+/* function for MMA's dual optimization of the approximant function */
+
+typedef struct {
+ int n; /* must be set on input to dimension of x */
+ const double *x, *lb, *ub, *sigma, *dfdx; /* arrays of length n */
+ const double *dfcdx; /* m-by-n array of fc gradients */
+ double fval, rho; /* must be set on input */
+ const double *fcval, *rhoc; /* arrays of length m */
+ double *xcur; /* array of length n, output each time */
+ double gval, wval, *gcval; /* output each time (array length m) */
+} dual_data;
+
+static double dual_func(int m, const double *y, double *grad, void *d_)
+{
+ dual_data *d = (dual_data *) d_;
+ int n = d->n;
+ const double *x = d->x, *lb = d->lb, *ub = d->ub, *sigma = d->sigma,
+ *dfdx = d->dfdx;
+ const double *dfcdx = d->dfcdx;
+ double rho = d->rho, fval = d->fval;
+ const double *rhoc = d->rhoc, *fcval = d->fcval;
+ double *xcur = d->xcur;
+ double *gcval = d->gcval;
+ int i, j;
+ double val;
+
+ val = d->gval = fval;
+ d->wval = 0;
+ for (i = 0; i < m; ++i) val += y[i] * (gcval[i] = fcval[i]);
+ if (grad) { for (i = 0; i < m; ++i) grad[i] = fcval[i]; }
+
+ for (j = 0; j < n; ++j) {
+ int x_pinned;
+ double a, b, dx, denom, denominv, ca, cb, sigma2;
+
+ /* first, compute xcur[j] for y */
+ a = dfdx[j];
+ b = fabs(dfdx[j]) * sigma[j] + 0.5 * rho;
+ for (i = 0; i < m; ++i) {
+ a += dfcdx[i*n + j] * y[i];
+ b += (fabs(dfcdx[i*n + j]) * sigma[j] + 0.5 * rhoc[i]) * y[i];
+ }
+ sigma2 = sigma[j] * sigma[j];
+ dx = -a * sigma2 / b;
+ xcur[j] = x[j] + dx;
+ if (xcur[j] > ub[j]) xcur[j] = ub[j];
+ if (xcur[j] < lb[j]) xcur[j] = lb[j];
+ if (xcur[j] > x[j] + 0.9*sigma[j]) xcur[j] = x[j] + 0.9*sigma[j];
+ if (xcur[j] < x[j] - 0.9*sigma[j]) xcur[j] = x[j] - 0.9*sigma[j];
+ x_pinned = xcur[j] != x[j] + dx;
+ dx = xcur[j] - x[j];
+
+ /* function value: */
+ denom = sigma2 - dx*dx;
+ denominv = 1.0 / denom;
+ ca = sigma2 * dx;
+ cb = dx*dx;
+ val += (a*ca + b*cb) * denominv;
+
+ /* update gval, wval, gcval */
+ d->gval += (dfdx[j] * ca + (fabs(dfdx[j])*sigma[j] + 0.5*rho) * cb)
+ * denominv;
+ d->wval += 0.5 * cb * denominv;
+ for (i = 0; i < m; ++i)
+ gcval[i] += (dfcdx[i*n+j] * ca
+ + (fabs(dfcdx[i*n+j])*sigma[j] + 0.5*rhoc[j]) * cb)
+ * denominv;
+
+ /* gradient */
+ if (grad)
+ for (i = 0; i < m; ++i) {
+ double dady, dbdy;
+ dady = dfcdx[i*n + j];
+ dbdy = fabs(dfcdx[i*n + j]) * sigma[j] + 0.5 * rhoc[i];
+ grad[i] += (dady*ca + dbdy*cb) * denominv;
+ if (!x_pinned) { /* dx/dy nonzero */
+ int k;
+ double dxdy;
+ dxdy = (a*dbdy - dady*b) * (sigma2 / (b*b));
+ grad[i] +=
+ ((sigma2 * dfdx[j] + 2 * (fabs(dfdx[j])*sigma[j]
+ + 0.5*rho) * dx)
+ * denom
+ + 2*dx * (dfdx[j]*ca + (fabs(dfdx[j])*sigma[j]
+ + 0.5*rho) * cb))
+ * (denominv*denominv * dxdy);
+ for (k = 0; k < m; ++k)
+ grad[i] +=
+ ((sigma2 * dfcdx[k*n + j]
+ + 2 * (fabs(dfcdx[k*n + j])*sigma[j]
+ + 0.5*rhoc[k]) * dx)
+ * denom
+ + 2*dx * (dfcdx[k*n + j]*ca
+ + (fabs(dfcdx[k*n + j])*sigma[j]
+ + 0.5*rhoc[k]) * cb))
+ * (denominv*denominv * dxdy) * y[k];
+ }
+ }
+ }
+
+ /* negate because we are maximizing */
+ if (grad) { for (i = 0; i < m; ++i) grad[i] = -grad[i]; }
+ return -val;
+}
+
+/***********************************************************************/
+