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