lib/minimize.c

   1 /* minimize - minimization for Bayes routines in XLISP-STAT and S      */
   2 /* XLISP-STAT 2.1 Copyright (c) 1990, by Luke Tierney                  */
   3 /* Additions to Xlisp 2.1, Copyright (c) 1989 by David Michael Betz    */
   4 /* You may give out copies of this software; for conditions see the    */
   5 /* file COPYING included with this distribution.                       */
   6
   7 /*
   8  * Nonlinear optimization modules adapted from Dennis and Schnabel,
   9  * "Numerical Methods for Unconstrained Optimization and Nonlinear
  10  * Equations."
  11  */
  12
  13 #include <stdio.h>
  14 #include "xmath.h"
  15 extern char buf[200];
  16 #define PRINTSTR(s) printf(s)
  17
  18 typedef struct {
  19   size_t n, k;
  20   void (*ffun)(), (*gfun)();
  21   double f, typf, new_f;
  22   double crit, new_crit;
  23   double *x, *new_x, *sx, *delf, *new_delf, *qnstep, *F;
  24   double **hessf, **H, **L, **new_delg;
  25   double gradtol, steptol, maxstep;
  26   int itncount, itnlimit, maxtaken, consecmax, retcode, termcode;
  27   int use_line_search, verbose, values_supplied, change_sign;
  28   double diagadd;
  29 } Iteration;
  30
  31
  32 extern double macheps();
  33
  34
  35 extern double choldecomp();
  36 static void eval_gradient(Iteration *);
  37 static void eval_next_gradient(Iteration *);
  38
  39
  40 /*
  41 This is a little more advanced section on functions, but is very useful. Take this for example:
  42   int applyeqn(int F(int), int max, int min) {
  43     int itmp;
  44
  45     itmp = F(int) + min;
  46     itmp = itmp - max;
  47     return itmp;
  48   }
  49
  50  What does this function do if we call it with applyeqn(square(x),
  51  y, z);? What happens is that the int F(int) is a reference to the
  52  function that is passed in as a parameter. Thus inside applyeqn where
  53  there is a call to F, it actually is a call to square! This is very
  54  useful if we have one set function, but wish to vary the input
  55  according to a particular function. So if we had a different function
  56  called cube we could change how we call applyeqn by calling the
  57  function by applyeqn(cube(x), y, z);.  */
  58
  59
  60 /************************************************************************/
  61 /**                                                                    **/
  62 /**                      Definitions and Globals                       **/
  63 /**                                                                    **/
  64 /************************************************************************/
  65
  66 # define nil 0L
  67 # define FALSE 0
  68 # define TRUE 1
  69
  70 # define INIT_GRAD_FRAC .001
  71 # define CONSEC_MAX_LIMIT 5
  72 # define ALPHA .0001
  73 # define MAX_STEP_FACTOR 1000
  74 # define GRADTOL_POWER 1.0 / 3.0
  75 # define STEPTOL_POWER 2.0 / 3.0
  76 # define ITNLIMIT 100
  77 # define VERBOSE 0
  78 # define USE_SEARCH TRUE
  79
  80 typedef double **RMatrix, *RVector;
  81
  82
  83 static char *termcodes[] = {"not yet terminated",
  84                             "gradient size is less than gradient tolerance",
  85                             "step size is less than step tolerance",
  86                             "no satisfactory step found in backtracking",
  87                             "iteration limit exceeded",
  88                             "maximum size step taken 5 iterations in a row"};
  89
  90 /************************************************************************/
  91 /**                                                                    **/
  92 /**                         Utility Functions                          **/
  93 /**                                                                    **/
  94 /************************************************************************/
  95
  96 static double Max(a, b)
  97         double a, b;
  98 {
  99   return(a > b ? a : b);
 100 }
 101
 102 static double Min(a, b)
 103         double a, b;
 104 {
 105   return(a > b ? b : a);
 106 }
 107
 108 /************************************************************************/
 109 /**                                                                    **/
 110 /**                    Cholesky Solving Functions                      **/
 111 /**                                                                    **/
 112 /************************************************************************/
 113
 114
 115 /* solve Ly = b for y */
 116 static void
 117 lsolve(size_t n, double *b, double **L, double *y)
 118 {
 119   size_t i, j;
 120
 121   for (i = 0; i < n; i++) {
 122     y[i] = b[i];
 123     for (j = 0; j < i; j++) y[i] -= L[i][j] * y[j];
 124     if (L[i][i] != 0) y[i] /= L[i][i];
 125   }
 126 }
 127
 128 /* solve (L^T)x = y for x */
 129 static void
 130 ltsolve(size_t n, double *y, double **L, double *x)
 131 {
 132   size_t i, j;
 133
 134   for (i = n - 1; i >= 0; i--) {
 135     x[i] = y[i];
 136     for (j = i + 1; j < n; j++) x[i] -= L[j][i] * x[j];
 137     if (L[i][i] != 0) x[i] /= L[i][i];
 138   }
 139 }
 140
 141 /* solve (L L^T) s = -g for s */
 142 static void
 143 cholsolve(size_t n,
 144           double *g, double **L,
 145           double *s)
 146 {
 147   size_t i;
 148
 149   /* solve Ly = g */
 150   lsolve(n, g, L, s);
 151
 152   /* solve L^Ts = y */
 153   ltsolve(n, s, L, s);
 154
 155   for (i = 0; i < n; i++) s[i] = -s[i];
 156 }
 157
 158 static void
 159 modelhess(size_t n, double *sx, double **H, double **L, double *diagadd)
 160 {
 161   size_t i, j;
 162   double sqrteps, maxdiag, mindiag, maxoff, maxoffl, maxposdiag, mu,
 163     maxadd, maxev, minev, offrow, sdd;
 164
 165   /* scale H on both sides with sx */
 166   for (i = 0; i < n; i++)
 167     for (j = 0; j < n; j++) H[i][j] /= sx[i] * sx[j];
 168
 169   /*
 170    * find mu to add to diagonal so no diagonal elements are negative
 171    * and largest diagonal dominates largest off diagonal element in H
 172    */
 173   sqrteps = sqrt(macheps());
 174   for (maxdiag = H[0][0], i = 1; i < n; i++)
 175     maxdiag = Max(maxdiag, H[i][i]);
 176   for (mindiag = H[0][0], i = 1; i < n; i++)
 177     mindiag = Min(mindiag, H[i][i]);
 178   maxposdiag = Max(0.0, maxdiag);
 179
 180   if (mindiag <= sqrteps * maxposdiag) {
 181     mu = 2 * (maxposdiag - mindiag) * sqrteps - mindiag;
 182     maxdiag += mu;
 183   }
 184   else mu = 0.0;
 185
 186   if (n > 1) {
 187     for (maxoff = fabs(H[0][1]), i = 0; i < n; i++)
 188       for (j = i + 1; j < n; j++)
 189         maxoff = Max(maxoff, fabs(H[i][j]));
 190   }
 191   else maxoff = 0.0;
 192
 193   if (maxoff * (1 + 2 * sqrteps) > maxdiag) {
 194     mu += (maxoff - maxdiag) + 2 * sqrteps * maxoff;
 195     maxdiag = maxoff * (1 + 2 * sqrteps);
 196   }
 197
 198   if (maxdiag == 0.0) {
 199     mu = 1;
 200     maxdiag = 1;
 201   }
 202
 203   if (mu > 0) for (i = 0; i < n; i++) H[i][i] += mu;
 204
 205   maxoffl = sqrt(Max(maxdiag, maxoff / n));
 206
 207   /*
 208    * compute the perturbed Cholesky decomposition
 209    */
 210   for (i = 0; i < n; i++)
 211     for (j = 0; j < n; j++) L[i][j] = H[i][j];
 212   choldecomp(L, n, maxoffl, &maxadd);
 213
 214   /*
 215    * add something to diagonal, if needed, to make positive definite
 216    * and recompute factorization
 217    */
 218   if (maxadd > 0) {
 219     maxev = H[0][0];
 220     minev = H[0][0];
 221     for (i = 0; i < n; i++) {
 222       for (offrow = 0.0, j = 0; j < n; j++)
 223         if (i != j) offrow += fabs(H[i][j]);
 224       maxev = Max(maxev, H[i][i] + offrow);
 225       minev = Min(minev, H[i][i] - offrow);
 226     }
 227     sdd = (maxev - minev) * sqrteps - minev;
 228     sdd = Max(sdd, 0.0);
 229     mu = Min(maxadd, sdd);
 230     for (i = 0; i < n; i++) H[i][i] += mu;
 231     for (i = 0; i < n; i++)
 232       for (j = 0; j < n; j++) L[i][j] = H[i][j];
 233     choldecomp(L, n, maxoffl, &maxadd);
 234     *diagadd = mu;
 235   }
 236   else *diagadd = 0.0;
 237
 238   /* unscale H and L */
 239   for (i = 0; i < n; i++)
 240     for (j = 0; j < n; j++) {
 241       H[i][j] *= sx[i] * sx[j];
 242       L[i][j] *= sx[i];
 243     }
 244 }
 245
 246 /************************************************************************/
 247 /**                                                                    **/
 248 /**                        Stopping Criteria                           **/
 249 /**                                                                    **/
 250 /************************************************************************/
 251
 252 static double
 253 gradsize(Iteration *iter, int new)
 254 {
 255   size_t n, i;
 256   double size, term, crit, typf;
 257   double *x, *delf, *sx;
 258
 259   n = iter->n + iter->k;
 260   crit = iter->crit;
 261   typf = iter->typf;
 262   x = iter->x;
 263   sx = iter->sx;
 264   delf = (new) ? iter->new_delf : iter->delf;
 265
 266   for (i = 0, size = 0.0; i < n; i++) {
 267     term = fabs(delf[i]) * Max(x[i], 1.0 / sx[i]) / Max(fabs(crit), typf);
 268     size = Max(size, term);
 269   }
 270   return(size);
 271 }
 272
 273 static double
 274 incrsize(Iteration *iter)
 275 {
 276   size_t n, i;
 277   double size, term;
 278   double *x, *new_x, *sx;
 279
 280   new_x = iter->new_x;
 281   n = iter->n + iter->k;
 282   x = iter->x;
 283   sx = iter->sx;
 284
 285   for (i = 0, size = 0.0; i < n; i++) {
 286     term = fabs(x[i] - new_x[i]) / Max(fabs(x[i]), 1.0 / sx[i]);
 287         size = Max(size, term);
 288   }
 289   return(size);
 290 }
 291
 292 static int /* AJR:TESTME guess! */
 293 stoptest0(Iteration *iter)
 294 {
 295   iter->consecmax = 0;
 296
 297   if (gradsize(iter, FALSE) <= INIT_GRAD_FRAC * iter->gradtol)
 298     iter->termcode = 1;
 299   else iter->termcode = 0;
 300
 301   return(iter->termcode);
 302 }
 303
 304 static int
 305 stoptest(Iteration *iter)
 306 {
 307   int termcode, retcode, itncount, itnlimit, maxtaken, consecmax;
 308   double gradtol, steptol;
 309
 310   retcode = iter->retcode;
 311   gradtol = iter->gradtol;
 312   steptol = iter->steptol;
 313   itncount = iter->itncount;
 314   itnlimit = iter->itnlimit;
 315   maxtaken = iter->maxtaken;
 316   consecmax = iter->consecmax;
 317
 318   termcode = 0;
 319   if (retcode == 1) termcode = 3;
 320   else if (gradsize(iter, TRUE) <= gradtol) termcode = 1;
 321   else if (incrsize(iter) <= steptol) termcode = 2;
 322   else if (itncount >= itnlimit) termcode = 4;
 323   else if (maxtaken) {
 324     consecmax++;
 325     if (consecmax >= CONSEC_MAX_LIMIT) termcode = 5;
 326   }
 327   else consecmax = 0;
 328
 329   iter->consecmax = consecmax;
 330   iter->termcode = termcode;
 331
 332   return(termcode);
 333 }
 334
 335 /************************************************************************/
 336 /**                                                                    **/
 337 /**               Function and Derivative Evaluation                   **/
 338 /**                                                                    **/
 339 /************************************************************************/
 340
 341 static void
 342 eval_funval(Iteration *iter)
 343 {
 344   int i;
 345
 346   (*(iter->ffun))(iter->x, &iter->f, nil, nil);
 347   if (iter->k == 0) iter->crit = iter->f;
 348   else {
 349     eval_gradient(iter);
 350     for (i = 0, iter->crit = 0.0; i < iter->n + iter->k; i++)
 351       iter->crit += 0.5 * pow(fabs(iter->delf[i]), 2.0);
 352   }
 353 }
 354
 355 static void
 356 eval_next_funval(Iteration *iter)
 357 {
 358   int i;
 359
 360   (*(iter->ffun))(iter->new_x, &iter->new_f, nil, nil);
 361   if (iter->k == 0) iter->new_crit = iter->new_f;
 362   else {
 363     eval_next_gradient(iter);
 364     for (i = 0, iter->new_crit = 0.0; i < iter->n + iter->k; i++)
 365       iter->new_crit += 0.5 * pow(fabs(iter->new_delf[i]), 2.0);
 366   }
 367 }
 368
 369 static void
 370 eval_gradient(Iteration *iter)
 371 {
 372   size_t i, j, n, k;
 373
 374   n = iter->n;
 375   k = iter->k;
 376
 377   (*(iter->ffun))(iter->x, nil, iter->delf, nil);
 378   if (k > 0) {
 379     (*(iter->gfun))(iter->x, iter->delf + n, nil, nil);
 380     (*(iter->gfun))(iter->x, nil, iter->new_delg, nil);
 381     for (i = 0; i < n; i++) {
 382       for (j = 0; j < k; j++)
 383         iter->delf[i] += iter->x[n + j] * iter->new_delg[j][i];
 384       for (j = 0; j < k; j++) {
 385         iter->hessf[n + j][i] = iter->new_delg[j][i];
 386         iter->hessf[i][n + j] = iter->new_delg[j][i];
 387       }
 388     }
 389   }
 390 }
 391
 392 static void
 393 eval_next_gradient(Iteration *iter)
 394 {
 395   size_t i, j, n, k;
 396
 397   n = iter->n;
 398   k = iter->k;
 399   (*(iter->ffun))(iter->new_x, nil, iter->new_delf, nil);
 400   if (k > 0) {
 401     (*(iter->gfun))(iter->new_x, iter->new_delf + n, nil, nil);
 402     (*(iter->gfun))(iter->new_x, nil, iter->new_delg, nil);
 403     for (i = 0; i < n; i++) {
 404       for (j = 0; j < k; j++)
 405         iter->new_delf[i] += iter->new_x[n + j] * iter->new_delg[j][i];
 406     }
 407   }
 408 }
 409
 410 static void eval_hessian(Iteration *iter)
 411 {
 412   size_t i, j, n, k;
 413
 414   n = iter->n;
 415   k = iter->k;
 416   (*(iter->ffun))(iter->x, nil, nil, iter->hessf);
 417   for (i = n; i < n + k; i++)
 418     for (j = n; j < n + k; j++) iter->hessf[i][j] = 0.0;
 419 }
 420
 421 /************************************************************************/
 422 /**                                                                    **/
 423 /**                     Backtracking Line Search                       **/
 424 /**                                                                    **/
 425 /************************************************************************/
 426
 427 static void
 428 linesearch(Iteration *iter)
 429 {
 430   size_t i, n;
 431   double newtlen, maxstep, initslope, rellength, lambda, minlambda,
 432     lambdatemp, lambdaprev, a, b, disc, critprev, f1, f2, a11, a12, a21, a22,
 433     del;
 434   double *qnstep, *delf, *x, *new_x, *sx;
 435
 436   n = iter->n + iter->k;
 437   if (! iter->use_line_search) {
 438     iter->maxtaken = FALSE;
 439     for (i = 0; i < n; i++)
 440       iter->new_x[i] = iter->x[i] + iter->qnstep[i];
 441     eval_next_funval(iter);
 442     iter->retcode = 0;
 443   }
 444   else{
 445     qnstep = iter->qnstep;
 446     maxstep = iter->maxstep;
 447     delf = (iter->k == 0) ? iter->delf : iter->F;
 448     x = iter->x;
 449     new_x = iter->new_x;
 450     sx = iter->sx;
 451
 452     iter->maxtaken = FALSE;
 453     iter->retcode = 2;
 454
 455     for (i = 0, newtlen = 0.0; i < n; i++)
 456       newtlen += pow(sx[i] * qnstep[i], 2.0);
 457     newtlen = sqrt(newtlen);
 458
 459     if (newtlen > maxstep) {
 460       for (i = 0; i < n; i++) qnstep[i] *= (maxstep / newtlen);
 461       newtlen = maxstep;
 462     }
 463
 464     for (i = 0, initslope = 0.0; i < n; i++) initslope += delf[i] * qnstep[i];
 465
 466     for (i = 0, rellength = 0.0; i < n; i++)
 467       rellength = Max(rellength, fabs(qnstep[i]) / Max(fabs(x[i]), 1.0 / sx[i]));
 468
 469     minlambda = (rellength == 0.0) ? 2.0 : iter->steptol / rellength;
 470
 471     lambda = 1.0;
 472         lambdaprev = 1.0; /* to make compilers happy */
 473         critprev = 1.0;   /* to make compilers happy */
 474     while (iter->retcode > 1) {
 475       for (i = 0; i < n; i++) new_x[i] = x[i] + lambda * qnstep[i];
 476       eval_next_funval(iter);
 477       if (iter->new_crit <= iter->crit + ALPHA * lambda * initslope) {
 478         iter->retcode = 0;
 479         if (lambda == 1.0 && newtlen > 0.99 * maxstep) iter->maxtaken = TRUE;
 480       }
 481       else if (lambda < minlambda) {
 482         iter->retcode = 1;
 483         iter->new_crit = iter->crit;
 484         for (i = 0; i < n; i++) new_x[i] = x[i];
 485       }
 486       else {
 487         if (lambda == 1.0) { /* first backtrack, quadratic fit */
 488           lambdatemp = - initslope
 489                      / (2 * (iter->new_crit - iter->crit - initslope));
 490         }
 491         else { /* all subsequent backtracks, cubic fit */
 492           del = lambda - lambdaprev;
 493           f1 = iter->new_crit - iter->crit - lambda * initslope;
 494           f2 = critprev - iter->crit - lambdaprev * initslope;
 495           a11 = 1.0 / (lambda * lambda);
 496           a12 = -1.0 / (lambdaprev * lambdaprev);
 497           a21 = -lambdaprev * a11;
 498           a22 = -lambda * a12;
 499           a = (a11 * f1 + a12 * f2) / del;
 500           b = (a21 * f1 + a22 * f2) / del;
 501           disc = b * b - 3 * a * initslope;
 502           if (a == 0) { /* cubic is a quadratic */
 503             lambdatemp = - initslope / (2 * b);
 504           }
 505           else { /* legitimate cubic */
 506             lambdatemp = (-b + sqrt(disc)) / (3 * a);
 507           }
 508           lambdatemp = Min(lambdatemp, 0.5 * lambda);
 509         }
 510         lambdaprev = lambda;
 511         critprev = iter->new_crit;
 512         lambda = Max(0.1 * lambda, lambdatemp);
 513         if (iter->verbose > 0) {
 514           sprintf(buf, "Backtracking: lambda = %g\n", lambda);
 515           PRINTSTR(buf);
 516         }
 517       }
 518     }
 519   }
 520 }
 521
 522 /************************************************************************/
 523 /**                                                                    **/
 524 /**                   Status Printing Functions                        **/
 525 /**                                                                    **/
 526 /************************************************************************/
 527
 528 static void
 529 print_header(Iteration *iter)
 530 {
 531   if (iter->verbose > 0) {
 532     sprintf(buf, "Iteration %d.\n", iter->itncount);
 533     PRINTSTR(buf);
 534   }
 535 }
 536
 537 static void
 538 print_status(Iteration *iter)
 539 {
 540   size_t i, j;
 541
 542   if (iter->verbose > 0) {
 543     sprintf(buf, "Criterion value = %g\n",
 544             (iter->change_sign) ? -iter->crit : iter->crit);
 545     PRINTSTR(buf);
 546     if (iter->verbose > 1) {
 547       PRINTSTR("Location = <");
 548       for (i = 0; i < iter->n + iter->k; i++) {
 549         sprintf(buf, (i < iter->n + iter->k - 1) ? "%g " : "%g>\n", iter->x[i]);
 550         PRINTSTR(buf);
 551       }
 552     }
 553     if (iter->verbose > 2) {
 554       PRINTSTR("Gradient = <");
 555       for (i = 0; i < iter->n + iter->k; i++) {
 556         sprintf(buf, (i < iter->n + iter->k - 1) ? "%g " : "%g>\n",
 557                 (iter->change_sign) ? -iter->delf[i] : iter->delf[i]);
 558         PRINTSTR(buf);
 559       }
 560     }
 561     if (iter->verbose > 3) {
 562       PRINTSTR("Hessian:\n");
 563       for (i = 0; i < iter->n + iter->k; i++) {
 564         for (j = 0; j < iter->n + iter->k; j++) {
 565           sprintf(buf, (j < iter->n + iter->k - 1) ? "%g " : "%g\n",
 566                   (iter->change_sign) ? -iter->hessf[i][j] : iter->hessf[i][j]);
 567           PRINTSTR(buf);
 568         }
 569       }
 570     }
 571   }
 572   if (iter->termcode != 0 && iter->verbose > 0) {
 573     sprintf(buf, "Reason for termination: %s.\n", termcodes[iter->termcode]);
 574     PRINTSTR(buf);
 575   }
 576 }
 577
 578 /************************************************************************/
 579 /**                                                                    **/
 580 /**                         Iteration Driver                           **/
 581 /**                                                                    **/
 582 /************************************************************************/
 583
 584 static void
 585 findqnstep(Iteration *iter)
 586 {
 587   size_t i, j, N, l;
 588
 589   if (iter->k == 0) {
 590     modelhess(iter->n, iter->sx, iter->hessf, iter->L, &iter->diagadd);
 591     cholsolve(iter->n, iter->delf, iter->L, iter->qnstep);
 592   }
 593   else {
 594     N = iter->n + iter->k;
 595     for (i = 0; i < N; i++) {
 596       for (l = 0, iter->F[i] = 0.0; l < N; l++)
 597         iter->F[i] += iter->hessf[i][l] * iter->delf[l];
 598       for (j = 0; j < N; j++)
 599         for (l = 0, iter->H[i][j] = 0.0; l < N; l++)
 600           iter->H[i][j] += iter->hessf[i][l] * iter->hessf[j][l];
 601     }
 602     modelhess(N, iter->sx, iter->H, iter->L, &iter->diagadd);
 603     cholsolve(N, iter->F, iter->L, iter->qnstep);
 604   }
 605 }
 606
 607 static void iterupdate(Iteration *iter)
 608 {
 609   size_t i, j, n, k;
 610
 611   n = iter->n;
 612   k = iter->k;
 613   iter->f = iter->new_f;
 614   iter->crit = iter->new_crit;
 615   for (i = 0; i < n + k; i++) {
 616     iter->x[i] = iter->new_x[i];
 617     iter->delf[i] = iter->new_delf[i];
 618   }
 619   for (i = 0; i < k; i++) {
 620     for (j = 0; j < n; j++) {
 621       iter->hessf[n + i][j] = iter->new_delg[i][j];
 622       iter->hessf[j][n + i] = iter->new_delg[i][j];
 623     }
 624   }
 625 }
 626
 627 static void
 628 mindriver(Iteration *iter)
 629 {
 630   iter->consecmax = 0;
 631   iter->itncount = 0;
 632   iter->termcode = 0;
 633   if (! iter->values_supplied) {
 634     eval_funval(iter);
 635     if (iter->k == 0) eval_gradient(iter);
 636     eval_hessian(iter);
 637   }
 638
 639   stoptest0(iter);
 640   print_header(iter);
 641   print_status(iter);
 642   while (iter->termcode == 0) {
 643     iter->itncount++;
 644     print_header(iter);
 645     findqnstep(iter);
 646     linesearch(iter);
 647     if (iter->k == 0) eval_next_gradient(iter);
 648     stoptest(iter);
 649     iterupdate(iter);
 650     eval_hessian(iter);
 651     print_status(iter);
 652   }
 653 }
 654
 655 /************************************************************************/
 656 /**                                                                    **/
 657 /**                   External Interface Routines                      **/
 658 /**                                                                    **/
 659 /************************************************************************/
 660
 661 static Iteration myiter;
 662
 663 size_t
 664 minworkspacesize(size_t n, size_t k)
 665 {
 666   size_t size;
 667
 668   /* x, new_x, sx, delf, new_delf, qnstep and F */
 669   size = 7 * sizeof(double) * (n + k);
 670
 671   /* hessf, H and L */
 672   size += 3 * (sizeof(double *) * (n + k)
 673                + sizeof(double) * (n + k) * (n + k));
 674
 675   /* delg and new_delg */
 676   size += 2 * (sizeof(double *) * k + sizeof(double) * n * k);
 677
 678   return(size);
 679 }
 680
 681 char *
 682 minresultstring(int code)
 683 {
 684   if (code <= 0) return("bad input data");
 685   else if (code <= 5) return(termcodes[code]);
 686   else return("unknown return code");
 687 }
 688
 689 void
 690 minsetup(size_t n, size_t k,
 691          void (*ffun)(), void (*gfun)(),
 692          double *x, double typf, double *typx, char *work)
 693      /*     int  (*ffun)(), (*gfun)();*/
 694 {
 695   Iteration *iter = &myiter;
 696   size_t i, j;
 697   double nx0, ntypx;
 698
 699   n = (n > 0) ? n : 0;
 700   k = (k > 0) ? k : 0;
 701
 702   iter->n = n;
 703   iter->k = k;
 704   iter->ffun = ffun; /* FIXME:AJR */
 705   iter->gfun = gfun;
 706
 707   iter->x = (double *) work; work += sizeof(double) * (n + k);
 708   iter->new_x = (double *) work; work += sizeof(double) * (n + k);
 709   iter->sx = (double *) work; work += sizeof(double) * (n + k);
 710   iter->delf = (double *) work; work += sizeof(double) * (n + k);
 711   iter->new_delf = (double *) work; work += sizeof(double) * (n + k);
 712   iter->qnstep = (double *) work; work += sizeof(double) * (n + k);
 713   iter->F = (double *) work; work += sizeof(double) * (n + k);
 714   for (i = 0; i < n; i++) {
 715     iter->x[i] = x[i];
 716     iter->sx[i] = (typx != nil && typx[i] > 0.0) ? 1.0 / typx[i] : 1.0;
 717   }
 718   for (i = 0; i < k; i++) {
 719     iter->x[n + i] = x[n + i];
 720     iter->sx[n + i] = 1.0;
 721   }
 722
 723   iter->hessf = (double **) work; work += sizeof(double *) * (n + k);
 724   for (i = 0; i < n + k; i++) {
 725     iter->hessf[i] = (double *) work;
 726     work += sizeof(double) * (n + k);
 727   }
 728   for (i = 0; i < n + k; i++)
 729     for (j = 0; j < n + k; j++) iter->hessf[i][j] = 0.0;
 730   iter->L = (double **) work; work += sizeof(double *) * (n + k);
 731   for (i = 0; i < n + k; i++) {
 732     iter->L[i] = (double *) work;
 733     work += sizeof(double) * (n + k);
 734   }
 735   iter->H = (double **) work; work += sizeof(double *) * (n + k);
 736   for (i = 0; i < n + k; i++) {
 737     iter->H[i] = (double *) work;
 738     work += sizeof(double) * (n + k);
 739   }
 740
 741   iter->new_delg = (k > 0) ? (double **) work : nil;
 742   work += sizeof(double *) * k;
 743   for (i = 0; i < k; i++) {
 744     iter->new_delg[i] = (double *) work;
 745     work += sizeof(double) * n;
 746   }
 747
 748   iter->typf = (typf > 0.0) ? typf : 1.0;
 749
 750   iter->verbose = VERBOSE;
 751   iter->use_line_search = USE_SEARCH;
 752   iter->itncount = 0;
 753   iter->itnlimit = ITNLIMIT;
 754   iter->gradtol = pow(macheps(), GRADTOL_POWER);
 755   iter->steptol = pow(macheps(), STEPTOL_POWER);
 756   for (i = 0, nx0 = 0.0, ntypx = 0.0; i < iter->n; i++) {
 757     nx0 += fabs(iter->new_x[i] / iter->sx[i]);
 758     ntypx += fabs(1.0 / iter->sx[i]);
 759   }
 760   iter->maxstep = MAX_STEP_FACTOR * Max(nx0, ntypx);
 761   iter->values_supplied = FALSE;
 762 }
 763
 764 void
 765 minsetoptions(double gradtol, double steptol, double maxstep,
 766               int itnlimit, int verbose, int use_search, int change_sign)
 767 {
 768   Iteration *iter = &myiter;
 769
 770   if (gradtol > 0.0) iter->gradtol = gradtol;
 771   if (steptol > 0.0) iter->steptol = steptol;
 772   if (maxstep > 0.0) iter->maxstep = maxstep;
 773   if (itnlimit >= 0) iter->itnlimit = itnlimit;
 774   if (verbose >= 0)  iter->verbose = verbose;
 775   iter->use_line_search = use_search;
 776   iter->change_sign = change_sign;
 777 }
 778
 779 void
 780 minsupplyvalues(double f, double *delf, double **hessf,
 781                 double *g, double **delg)
 782 {
 783   Iteration *iter = &myiter;
 784   size_t i, j, n, k;
 785
 786   n = iter->n;
 787   k = iter->k;
 788
 789   iter->f = f;
 790   for (i = 0; i < n; i++) {
 791     iter->delf[i] = delf[i];
 792     for (j = 0; j < k; j++)
 793       iter->delf[i] += iter->x[n + j] * delg[j][i];
 794     for (j = 0; j < n; j++) iter->hessf[i][j] = hessf[i][j];
 795   }
 796   for (i = 0; i < k; i++) {
 797     iter->delf[n + i] = g[i];
 798     for (j = 0; j < n; j++) {
 799       iter->hessf[n + i][j] = delg[i][j];
 800       iter->hessf[j][n + i] = delg[i][j];
 801     }
 802   }
 803   if (k == 0) iter->crit = f;
 804   else {
 805     for (i = 0, iter->crit = 0.0; i < n + k; i++)
 806       iter->crit += 0.5 * pow(fabs(iter->delf[i]), 2.0);
 807   }
 808   iter->values_supplied = TRUE;
 809 }
 810
 811 void
 812 minimize()
 813 {
 814   mindriver(&myiter);
 815 }
 816
 817 void
 818 minresults(double *x, double *pf, double *pcrit, double *delf, double **hessf,
 819            double *g, double **delg, int *pcount, int *ptermcode, double *diagadd)
 820 {
 821   Iteration *iter = &myiter;
 822   size_t i, j, n, k;
 823
 824   n = iter->n;
 825   k = iter->k;
 826
 827   if (pf != nil) *pf = iter->f;
 828   if (pcrit != nil) *pcrit = iter->crit;
 829   for (i = 0; i < n; i++) {
 830     if (x != nil) x[i] = iter->x[i];
 831     if (delf != nil) delf[i] = iter->delf[i];
 832     for (j = 0; j < n; j++) if (hessf != nil) hessf[i][j] = iter->hessf[i][j];
 833   }
 834   for (i = 0; i < k; i++) {
 835     if (x != nil) x[n + i] = iter->x[n + i];
 836     if (g != nil) g[i] = iter->delf[n + i];
 837     for (j = 0; j < n; j++)
 838       if (delg != nil) delg[i][j] = iter->hessf[n + i][j];
 839   }
 840   if (pcount != nil) *pcount = iter->itncount;
 841   if (ptermcode != nil) *ptermcode = iter->termcode;
 842   if (diagadd != nil) *diagadd = iter->diagadd;
 843 }
 844
 845
 846 double pdlogdet(int n, double **H)
 847 {
 848   int i;
 849   double logdet, maxadd;
 850
 851   choldecomp(H, n, 0.0, &maxadd);
 852   for (i = 0, logdet = 0.0; i < n; i++) logdet += 2.0 * log(H[i][i]);
 853   return(logdet);
 854 }
 855
 856
 857 /*
 858   TODO
 859
 860   return amount added to make pos definite
 861
 862   scaling for constraints
 863
 864   alternate global strategies
 865
 866   callback function for verbose mode
 867 */
 868