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