lib/functions.c

   1 /* functions - callbacks 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 #include <stdio.h>
   8 #include "xmath.h"
   9
  10 #define PRINTSTR(s) printf(s)
  11
  12 extern void *S_alloc(), *calloc(), *realloc();
  13 extern double macheps();
  14 char *minresultstring();
  15
  16 extern void choldecomp();
  17
  18 extern void minsupplyvalues(double , double *, double **,
  19                             double *, double **);
  20 extern void minimize();
  21 extern void minresults(double *, double *, double *, double *, double **,
  22                        double *, double **, int *, int *, double *);
  23 extern void minsetup(int n, int k,
  24                      int *ffun(), int *gfun(),
  25                      double *x, double typf, double *typx, char *work);
  26                          /*     int  (*ffun)(), (*gfun)();*/
  27 extern void minsetoptions(double , double, double, int , int , int , int);
  28
  29 extern void choldecomp();
  30
  31
  32
  33 /* next 2 from cbayes.c */
  34 extern void Recover(char *, char *);
  35 extern void call_S(char *fun, long narg, char **args, char **mode, long *length,char **names,
  36                    long nvals, char **values);
  37
  38 /* next 2 from derivatives.c */
  39 /*
  40 extern void numergrad(int n, double *x, double *grad, double *fsum,
  41                       int *ffun(), double h, double *typx);
  42 extern void numerhess(int n, double *x, double **hess, double f,
  43                       double *fsum, void ffun(),
  44                       double h, double *typx);
  45 */
  46 extern void numergrad(int , double *, double *, double *,
  47                       int ffun(double *, double *, double *, double **),
  48                       double , double *);
  49 extern void numerhess(int , double *, double **, double ,
  50                       double *,
  51                       int ffun(double *, double *, double *, double **),
  52                       double , double *);
  53
  54 /* next from minimize */
  55 extern int minworkspacesize(int, int);
  56
  57 /************************************************************************/
  58 /**                                                                    **/
  59 /**                      Definitions and Globals                       **/
  60 /**                                                                    **/
  61 /************************************************************************/
  62
  63 /* #define NULL 0L  already defined in stddef */
  64
  65 #define nil 0L
  66 #define TRUE 1
  67 #define FALSE 0
  68
  69 #define ROOT2PI 2.50662827463100050241
  70 #define PI_INV .31830988618379067153
  71
  72 #define GRADTOL_POWER 1.0 / 3.0
  73 #define H_POWER 1.0 / 6.0
  74
  75 /*typedef double **RMatrix, *RVector;*/
  76
  77 typedef struct{
  78   char *f, **sf, **g;
  79   int n, k;
  80   int change_sign, fderivs;
  81   int *gderivs;
  82   double typf, h, dflt;
  83   double *typx, *fsum, *cvals, *ctarget;
  84   double **gfsum;
  85 } Fundata;
  86
  87 static Fundata func, gfuncs, cfuncs;
  88
  89 /************************************************************************/
  90 /**                                                                    **/
  91 /**                         Memory Utilities                           **/
  92 /**                                                                    **/
  93 /************************************************************************/
  94
  95 /* this function is used to maintain a statically allocated piece of    */
  96 /* memory of a specified size. If a larger piece is needed the pointer  */
  97 /* is realloced. This allows functions using memory allocation to be    */
  98 /* called repeatedly (but not recursively) from within the same call    */
  99 /* from S. It attempts to avoid the danger of dangling callocs.         */
 100
 101 void static
 102 makespace(pptr, size)
 103      char **pptr;
 104      int size;
 105 {
 106   if (size <= 0) return;
 107   if (*pptr == nil) *pptr = calloc(size, 1);
 108   else *pptr = realloc(*pptr, size);
 109   if (size > 0 && *pptr == nil) Recover("memory allocation failed", NULL);
 110 }
 111
 112 /************************************************************************/
 113 /**                                                                    **/
 114 /**                    Functions Evaluation Routines                   **/
 115 /**                                                                    **/
 116 /************************************************************************/
 117
 118 /*
 119  * All Hessian evaluations by numerical derivatives assume the gradient is
 120  * evaluated first at the same location. The results are cached away.
 121  */
 122
 123 /* install log posterior function */
 124 static void
 125 install_func(char *f, char **sf, int n, int change_sign, /*?? */
 126                     double typf, double h, double *typx, double dflt)
 127 {
 128   int i;
 129   static int inited = FALSE;
 130
 131   if (! inited) {
 132     func.typx = nil;
 133     func.fsum = nil;
 134     inited = TRUE;
 135   }
 136   makespace(&func.typx, n * sizeof(double));
 137   makespace(&func.fsum, n * sizeof(double));
 138
 139   func.f = f;
 140   func.sf = sf;
 141   func.n = n;
 142   func.change_sign = change_sign;
 143   func.typf = (typf > 0.0) ? typf : 1.0;
 144   func.h = (h > 0.0) ? h : pow(macheps(), H_POWER);
 145   for (i = 0; i < n; i++)
 146     func.typx[i] = (typx != nil && typx[i] > 0.0) ? typx[i] : 1.0;
 147   func.dflt = dflt;
 148   func.fderivs = 0;
 149 }
 150
 151 /* install tilt functions */
 152 static void
 153 install_gfuncs(char **g, int n, int k, int change_sign, double h, double *typx)
 154 {
 155   int i;
 156   static int inited = FALSE;
 157   static double *gfsumdata = nil;
 158
 159   if (! inited) {
 160     gfuncs.typx = nil;
 161     gfuncs.gfsum = nil;
 162     gfuncs.gderivs = nil;
 163     inited = TRUE;
 164   }
 165   makespace(&gfuncs.typx, n * sizeof(double));
 166   makespace(&gfuncs.gfsum, k * sizeof(double *));
 167   makespace(&gfsumdata, k * n * sizeof(double));
 168   makespace(&gfuncs.gderivs, k *sizeof(int));
 169
 170   gfuncs.g = g;
 171   gfuncs.n = n;
 172   gfuncs.k = k;
 173   gfuncs.change_sign = change_sign;
 174   gfuncs.h = (h > 0.0) ? h : pow(macheps(), H_POWER);
 175   for (i = 0; i < n; i++)
 176     gfuncs.typx[i] = (typx != nil && typx[i] > 0.0) ? typx[i] : 1.0;
 177   for (i = 0; i < k; i++) gfuncs.gfsum[i] = gfsumdata + i * n;
 178   return;
 179 }
 180
 181 /* install constraint functions */
 182 static void
 183 install_cfuncs(char **g, int n, int k, double *ctarget, double h, double *typx)
 184 {
 185   int i;
 186   static int inited = FALSE;
 187
 188   if (! inited) {
 189     cfuncs.typx = nil;
 190     cfuncs.fsum = nil;
 191     cfuncs.gderivs = nil;
 192     inited = TRUE;
 193   }
 194   makespace(&cfuncs.typx, n * sizeof(double));
 195   makespace(&cfuncs.fsum, n * sizeof(double));
 196   makespace(&cfuncs.gderivs, k * sizeof(int));
 197
 198   cfuncs.g = g;
 199   cfuncs.n = n;
 200   cfuncs.k = k;
 201   cfuncs.h = (h > 0.0) ? h : pow(macheps(), H_POWER);
 202   for (i = 0; i < n; i++)
 203     cfuncs.typx[i] = (typx != nil && typx[i] > 0.0) ? typx[i] : 1.0;
 204   cfuncs.ctarget = ctarget;
 205   return;
 206 }
 207
 208 static int
 209 in_support(char **ff, int n, double *x)
 210 {
 211   char *args[1], *values[1];
 212   int *result;
 213   char *mode[1];
 214   long length[1];
 215
 216   if (ff == nil || ff[0] == nil) {
 217     return(TRUE);
 218   } else {
 219     mode[0] = "double";
 220     length[0] =n;
 221     args[0] = (char *) x;
 222     call_S(ff[0], 1L, args, mode, length, 0L, 1L, values);
 223     result = (int *) values[0];
 224     return(result[0]);
 225   }
 226 }
 227
 228 /* callback for logposterior evaluation */
 229 static int
 230 evalfunc(double *x, double *pval, double *grad, double **hess)
 231 {
 232   char *args[1], *values[3];
 233   double *result, val;
 234   char *mode[1];
 235   long length[1];
 236   int i, j;
 237
 238   for (i = 0; i < 3; i++) values[i] = nil;
 239
 240   if (in_support(func.sf, func.n, x)) {
 241     if (pval != nil || func.fderivs > 0 || hess != nil) {
 242       mode[0] = "double";
 243       length[0] = func.n;
 244       args[0] = (char *) x;
 245       call_S(func.f, 1L, args, mode, length, 0L, 3L, values);
 246       result = (double *) values[0];
 247           val = (! func.change_sign) ? result[0] : -result[0];
 248       if (pval != nil) *pval = val;
 249       if (values[2] != nil) func.fderivs = 2;
 250       else if (values[1] != nil) func.fderivs = 1;
 251       else func.fderivs = 0;
 252     }
 253     if (grad != nil) {
 254       if (func.fderivs > 0) {
 255         result = (double *) values[1];
 256         for (i = 0; i < func.n; i++)
 257           grad[i] = (! func.change_sign) ? result[i] : -result[i];
 258       }
 259       else {
 260         numergrad(func.n, x, grad, func.fsum, evalfunc, func.h, func.typx);
 261       }
 262     }
 263     if (hess != nil) {
 264       if (func.fderivs == 2) {
 265         result = (double *) values[2];
 266         for (i = 0; i < func.n; i++)
 267           for (j = 0; j < func.n; j++)
 268             hess[i][j] = (! func.change_sign) ? result[i + j * func.n]
 269                                               : -result[i + j * func.n];
 270       }
 271       else {
 272         if (func.fderivs == 1) /* kludge to get fsum for analytic gradients */
 273           numergrad(func.n, x, func.fsum, func.fsum,
 274                     evalfunc, func.h, func.typx);
 275         numerhess(func.n, x, hess, val, func.fsum, evalfunc, func.h, func.typx);
 276       }
 277     }
 278     return(TRUE);
 279   } else {
 280     if (pval != nil) {
 281       *pval = func.dflt;
 282     }
 283     return(FALSE);
 284   }
 285   return(TRUE);
 286 }
 287
 288
 289 /* callback for tilt function evaluation */
 290 static int which_gfunc;
 291
 292 static int
 293 evalgfunc(double *x, double *pval, double *grad, double **hess)
 294 {
 295   char *args[1], *values[3];
 296   double *result, val;
 297   char *mode[1];
 298   long length[1];
 299   int i, j;
 300
 301   for (i = 0; i < 3; i++) values[i] = nil;
 302
 303   if (pval != nil || gfuncs.gderivs[which_gfunc] > 0 || hess != nil) {
 304     mode[0] = "double";
 305     length[0] = gfuncs.n;
 306     args[0] = (char *) x;
 307     call_S(gfuncs.g[which_gfunc], 1L, args, mode, length, 0L, 3L, values);
 308     result = (double *) values[0];
 309         val = result[0];
 310     if (pval != nil) *pval = result[0];
 311     if (values[2] != nil) gfuncs.gderivs[which_gfunc] = 2;
 312     else if (values[1] != nil) gfuncs.gderivs[which_gfunc] = 1;
 313     else gfuncs.gderivs[which_gfunc] = 0;
 314   }
 315   if (grad != nil) {
 316     if (gfuncs.gderivs[which_gfunc] > 0) {
 317       result = (double *) values[1];
 318       for (i = 0; i < gfuncs.n; i++) grad[i] = result[i];
 319     }
 320     else {
 321       numergrad(gfuncs.n, x, grad, gfuncs.gfsum[which_gfunc], evalgfunc,
 322                 gfuncs.h, gfuncs.typx);
 323     }
 324   }
 325   if (hess != nil) {
 326     if (gfuncs.gderivs[which_gfunc] == 2) {
 327       result = (double *) values[2];
 328       for (i = 0; i < gfuncs.n; i++)
 329         for (j = 0; j < gfuncs.n; j++)
 330           hess[i][j] = result[i + j * gfuncs.n];
 331     }
 332     else {
 333       /* kludge to get fsum if analytic gradient used */
 334       if (gfuncs.gderivs[which_gfunc] == 1)
 335         numergrad(gfuncs.n, x, gfuncs.gfsum[which_gfunc],
 336                   gfuncs.gfsum[which_gfunc], evalgfunc, gfuncs.h, gfuncs.typx);
 337       numerhess(gfuncs.n, x, hess, val, gfuncs.gfsum[which_gfunc], evalgfunc,
 338                 gfuncs.h, gfuncs.typx);
 339     }
 340   }
 341   return(TRUE);
 342 }
 343
 344 /* callback for constraint function evaluation */
 345 static int which_cfunc;
 346
 347 static int
 348 evalcfunc(double *x, double *pval, double *grad, double **hess)
 349 {
 350   char *args[1], *values[3];
 351   double *result, val;
 352   char *mode[1];
 353   long length[1];
 354   int i, j;
 355
 356   if (pval != nil || cfuncs.gderivs[which_cfunc] > 0 || hess != nil) {
 357     mode[0] = "double";
 358     length[0] = cfuncs.n;
 359     args[0] = (char *) x;
 360     call_S(cfuncs.g[which_cfunc], 1L, args, mode, length, 0L, 3L, values);
 361     result = (double *) values[0];
 362         val = result[0];
 363     if (pval != nil) {
 364       *pval = result[0];
 365       if (cfuncs.ctarget != nil) *pval -= cfuncs.ctarget[which_cfunc];
 366     }
 367     if (values[2] != nil) cfuncs.gderivs[which_cfunc] = 2;
 368     else if (values[1] != nil) cfuncs.gderivs[which_cfunc] = 1;
 369     else cfuncs.gderivs[which_cfunc] = 0;
 370   }
 371   if (grad != nil) {
 372     if (cfuncs.gderivs[which_cfunc] > 0) {
 373       result = (double *) values[1];
 374       for (i = 0; i <cfuncs.n; i++) grad[i] = result[i];
 375     }
 376     else {
 377       numergrad(cfuncs.n, x, grad, cfuncs.fsum, evalcfunc,
 378                 cfuncs.h, cfuncs.typx);
 379     }
 380   }
 381   if (hess != nil) {
 382     if (cfuncs.gderivs[which_cfunc] == 2) {
 383       result = (double *) values[2];
 384       for (i = 0; i <cfuncs.n; i++)
 385         for (j = 0; j <cfuncs.n; j++)
 386           hess[i][j] = result[i + j * cfuncs.n];
 387     }
 388     else {
 389       /* kludge to get fsum if analytic gradient used */
 390       if (cfuncs.gderivs[which_cfunc] == 1)
 391         numergrad(cfuncs.n, x, cfuncs.fsum, cfuncs.fsum, evalcfunc,
 392                   cfuncs.h, cfuncs.typx);
 393       numerhess(cfuncs.n, x, hess, val, cfuncs.fsum, evalcfunc,
 394                 cfuncs.h, cfuncs.typx);
 395     }
 396   }
 397   return(TRUE);
 398 }
 399
 400 /* S front end for logposterior evaluation */
 401 void
 402 evalfront(char **ff, int *n, double *x, double *val, double *grad,
 403           double *phess, double *h, double *typx)
 404 {
 405   int i;
 406   static double **hess = nil;
 407
 408   install_func(ff[0], nil, *n, FALSE, 1.0, *h, typx, 0.0);
 409   if (phess == nil) hess = nil;
 410   else {
 411     makespace(&hess, *n * sizeof(double *));
 412     for (i = 0; i < *n; i++, phess += *n) hess[i] = phess;
 413   }
 414   evalfunc(x, val, grad, hess);
 415 }
 416
 417 /* S front end for tilt function evaluation */
 418 void
 419 gevalfront(char **gg, int *n, int *m, double *x, double *h,
 420            double *typx, double *val, double *grad)
 421 {
 422   int i;
 423
 424   install_gfuncs(gg, *n, *m, FALSE, *h, typx);
 425   for (i = 0; i < *m; i++, val++) {
 426     which_gfunc = i;
 427     evalgfunc(x, val, grad, nil);
 428     if (grad != nil) grad += *n;
 429   }
 430   return;
 431 }
 432
 433 /************************************************************************/
 434 /**                                                                    **/
 435 /**                     Derivative Scaling Routines                    **/
 436 /**                                                                    **/
 437 /************************************************************************/
 438
 439 /*
 440 void
 441 derivscalefront(char **ff, int *n, double *x, double *h, double *typx, double *tol, int *info)
 442 {
 443   int i;
 444
 445   if (*tol <= 0.0) *tol = pow(macheps(), GRADTOL_POWER);
 446
 447   install_func(ff[0], nil, *n, TRUE, 1.0, *h, typx, 0.0);
 448   *info = check_derivs(x, *tol);
 449
 450   *h = func.h;
 451   for (i = 0; i < *n; i++) typx[i] = func.typx[i];
 452   return;
 453 }
 454
 455
 456 static int
 457 check_derivs(double *x, double drvtol)
 458 {
 459   static double *grad = nil, work = nil;
 460   static double **hess = nil;
 461   int i, error;
 462
 463   grad = (double *) S_alloc(func.n, sizeof(double));
 464   hess = (double **) S_alloc(func.n, sizeof(double *));
 465   work = (double *) S_alloc(func.n + func.n * func.n, sizeof(double));
 466
 467   for (i = 0; i < func.n; i++) {
 468     hess[i] = work;
 469     work += func.n;
 470   }
 471
 472   error = derivscale(func.n, x, grad, hess, func.fsum, evalfunc,
 473                      &func.h, func.typx, drvtol, work);
 474   return(error);
 475 }
 476 */
 477
 478 /************************************************************************/
 479 /**                                                                    **/
 480 /**                    Importance Sampling Routines                    **/
 481 /**                                                                    **/
 482 /************************************************************************/
 483
 484 /* joint density of normal-cauchy mixture */
 485 static double
 486 dncmix(double *x, int n, double p)
 487 {
 488   int i;
 489   double dens;
 490
 491   for (i = 0, dens = 1.0; i < n; i++) {
 492     dens *= p * exp(-0.5 * x[i] * x[i]) / ROOT2PI
 493       + (1 - p) * PI_INV / (1.0 + x[i] * x[i]);
 494   }
 495   return(dens);
 496 }
 497
 498 /*
 499  * S front end for computing sample from transformed normal-cauchy
 500  * mixture and importance sampling weights
 501  */
 502 void
 503 samplefront(char **ff, char **sf, char **rf,
 504             double *p, int *n,
 505             double *x, double *ch, int *N, double *y, double *w)
 506 {
 507   double val;
 508   int i, j, k;
 509   char *args[2], *values[1];
 510   double *result, mval, c, dens, m;
 511   char *mode[2];
 512   long length[2];
 513
 514   /* get the random variables */
 515   mode[0] = "double"; mode[1] = "double";
 516   length[0] = 1; length[1] = 1;
 517   m = *N * *n; args[0] = (char *) &m; args[1] = (char *) p;
 518   call_S(rf[0], 2L, args, mode, length, 0L, 1L, values);
 519   result = (double *) values[0];
 520   for (i = 0; i < m; i++) y[i] = result[i];
 521
 522   /* construct the sample and the weights */
 523   install_func(ff[0], sf, *n, FALSE, 1.0, -1.0, nil, 0.0);
 524   c = 1.0 / pow(ROOT2PI, (double) *n);
 525   evalfunc(x, &mval, nil, nil);
 526   for (i = 0; i < *N; i++, y += *n) {
 527     dens = dncmix(y, *n, *p);
 528     for (j = 0; j < *n; j++) {
 529       val = x[j];
 530       for (k = j; k < *n; k++) val += y[k] * ch[j + *n * k];
 531       y[j] = val;
 532     }
 533     if (evalfunc(y, &val, nil, nil)) w[i] = exp(val - mval) * c / dens;
 534     else w[i] = 0.0;
 535   }
 536   return;
 537 }
 538
 539
 540 /************************************************************************/
 541 /**                                                                    **/
 542 /**                       Maximization Routines                        **/
 543 /**                                                                    **/
 544 /************************************************************************/
 545
 546 typedef struct {
 547   int n, m, k, itnlimit, backtrack, verbose, vals_suppl, exptilt;
 548   int count, termcode;
 549 } MaxIPars;
 550
 551 typedef struct {
 552   double typf, h, gradtol, steptol, maxstep, dflt, tilt, newtilt, hessadd;
 553 } MaxDPars;
 554
 555 struct {
 556   double tilt;
 557   double *gval;
 558   double  **ggrad, **ghess;
 559   int exptilt;
 560   double *tscale;
 561 } tiltinfo;
 562
 563 static void
 564 add_tilt(double *x, double *pval, double *grad, double **hess,
 565          double tilt, int exptilt)
 566 {
 567   int i, j, k, n = func.n, m = gfuncs.k;
 568   double *gval, *ggrad, **ghess, etilt;
 569
 570   if (m == 0) return;
 571
 572   if (gfuncs.change_sign) tilt = -tilt;
 573
 574   for (k = 0; k < m; k++) {
 575     gval = (pval != nil) ? tiltinfo.gval + k : nil;
 576     ggrad = (grad != nil) ? tiltinfo.ggrad[k] : nil;
 577     ghess = (hess != nil) ? tiltinfo.ghess : nil;
 578
 579     which_gfunc = k;
 580     evalgfunc(x, gval, ggrad, ghess);
 581
 582     if (exptilt) {
 583       etilt = (tiltinfo.tscale != nil) ? tilt / tiltinfo.tscale[k] : tilt;
 584       if (pval != nil) *pval += etilt * *gval;
 585       if (grad != nil)
 586         for (i = 0; i < n; i++) grad[i] += etilt * ggrad[i];
 587       if (hess != nil)
 588         for (i = 0; i < n; i++)
 589           for (j = 0; j < n; j++) hess[i][j] += etilt * ghess[i][j];
 590     }
 591     else {
 592       gval = tiltinfo.gval;
 593       ggrad = tiltinfo.ggrad[k];
 594       ghess = tiltinfo.ghess;
 595       if (gval[k] <= 0.0) Recover("nonpositive function value", NULL);
 596       if (pval != nil) *pval += tilt * log(gval[k]);
 597       if (grad != nil)
 598         for (i = 0; i < n; i++) grad[i] += tilt * ggrad[i] / gval[k];
 599       if (hess != nil)
 600         for (i = 0; i < n; i++)
 601           for (j = 0; j < n; j++)
 602             hess[i][j] +=
 603               tilt * (ghess[i][j] / gval[k]
 604                       - (ggrad[i] / gval[k]) * (ggrad[j] / gval[k]));
 605     }
 606   }
 607 }
 608
 609 static void
 610 set_tilt_info(int n, int m,
 611               double tilt, int exptilt, double *tscale)
 612 {
 613   static double *hessdata = nil, *graddata = nil;
 614   int i;
 615   static int inited = FALSE;
 616
 617   if (! inited) {
 618     tiltinfo.gval = nil;
 619     tiltinfo.ggrad = nil;
 620     tiltinfo.ghess = nil;
 621     inited = TRUE;
 622   }
 623   makespace(&tiltinfo.gval, n * sizeof(double));
 624   makespace(&tiltinfo.ggrad, m * sizeof(double *));
 625   makespace(&tiltinfo.ghess, n * sizeof(double *));
 626   makespace(&graddata, n * m * sizeof(double));
 627   makespace(&hessdata, n * n * sizeof(double));
 628
 629   tiltinfo.tilt = tilt;
 630   tiltinfo.exptilt = exptilt;
 631   for (i = 0; i < m; i++) tiltinfo.ggrad[i] = graddata + i * n;
 632   for (i = 0; i < n; i++) tiltinfo.ghess[i] = hessdata + i * n;
 633   tiltinfo.tscale = tscale;
 634 }
 635
 636
 637 static void
 638 minfunc(double *x, double *pval, double *grad, double **hess)
 639 {
 640   int k = gfuncs.k;
 641
 642   if (evalfunc(x, pval, grad, hess) && (k > 0))
 643     add_tilt(x, pval, grad, hess, tiltinfo.tilt, tiltinfo.exptilt);
 644 }
 645
 646 void
 647 constfunc(double *x, double *vals, double **jac, double **hess)
 648 {
 649   int i, k = cfuncs.k;
 650   double *pvali, *jaci;
 651
 652   for (i = 0; i < k; i++) {
 653     pvali = (vals != nil) ? vals + i : nil;
 654     jaci = (jac != nil) ? jac[i] : nil;
 655     which_cfunc = i;
 656     evalcfunc(x, pvali, jaci, nil);
 657   }
 658 }
 659
 660 void
 661 maxfront(char **ff, char **gf, char **cf,
 662          double *x, double *typx, double *fvals, double *gvals, double *cvals, double *ctarget,
 663          MaxIPars *ipars, MaxDPars *dpars,
 664          double *tscale, char **msg)
 665 {
 666   static char *work = nil;
 667   static double **H = nil, **cJ = nil;
 668   double *pf, *grad, *c;
 669   int i, n, m, k;
 670   void (*cfun)();
 671
 672   if (ipars->verbose > 0) PRINTSTR("maximizing...\n");
 673
 674   n = ipars->n;
 675   m = ipars->m;
 676   k = ipars->k;
 677   if (k >= n) Recover("too many constraints", NULL);
 678
 679   makespace(&H, n * sizeof(double *));
 680   makespace(&work, minworkspacesize(n, k));
 681
 682   pf = fvals; fvals++;
 683   grad = fvals; fvals += n;
 684   for (i = 0; i < n; i++, fvals += n) H[i] = fvals;
 685   set_tilt_info(n, m, dpars->newtilt, ipars->exptilt, tscale);
 686
 687   if (k == 0) {
 688     c = nil;
 689     cJ = nil;
 690     cfun = nil;
 691   } else {
 692     c = cvals;
 693     cvals += k;
 694     makespace(&cJ, k * sizeof(double *));
 695     for (i = 0; i < k; i++) cJ[i] = cvals + i * n;
 696     cfun = &constfunc; /* AJR: pointer to constfunc? */
 697   }
 698
 699   install_func(ff[0], nil, n, TRUE, dpars->typf, dpars->h, typx, dpars->dflt);
 700   install_gfuncs(gf, n, m, TRUE, dpars->h, typx);
 701   install_cfuncs(cf, n, k, ctarget, dpars->h, typx);
 702
 703   minsetup(n, k, &minfunc, &cfun, x, dpars->typf, typx, work); /* AJR: FIXME arg 3,4 */
 704   minsetoptions(dpars->gradtol, dpars->steptol, dpars->maxstep,
 705                 ipars->itnlimit, ipars->verbose, ipars->backtrack, TRUE);
 706
 707   if (ipars->vals_suppl) {
 708     for (i = 0; i < k; i++) c[i] -= ctarget[i];
 709     if (dpars->newtilt != dpars->tilt) {
 710       add_tilt(x, pf, grad, H, dpars->newtilt - dpars->tilt, ipars->exptilt);
 711       dpars->tilt = dpars->newtilt;
 712     }
 713     minsupplyvalues(*pf, grad, H, c, cJ);
 714   }
 715
 716   minimize();
 717   minresults(x, pf, nil, grad, H, c, cJ, &ipars->count, &ipars->termcode,
 718              &dpars->hessadd);
 719   msg[0] = minresultstring(ipars->termcode);
 720
 721   for (i = 0; i < k; i++) c[i] += ctarget[i];
 722   ipars->vals_suppl = TRUE;
 723 }
 724
 725 /************************************************************************/
 726 /**                                                                    **/
 727 /**                     Log Laplace Approximation                      **/
 728 /**                                                                    **/
 729 /************************************************************************/
 730
 731 void
 732 loglapdet(double *fvals, double *cvals, MaxIPars *ipars, MaxDPars *dpars,
 733           double *val, int *detonly)
 734 {
 735   int i, j, l, n = ipars->n, k = ipars->k;
 736   double f = -fvals[0], *hessdata = fvals + n + 1, *cgraddata = cvals + k;
 737   double ldL, ldcv, maxadd;
 738   static double **hess = nil, **cgrad = nil;
 739
 740   if (k >= n) Recover("too many constraints", NULL);
 741
 742   makespace(&hess, n * sizeof(double *));
 743   makespace(&cgrad, k * sizeof(double *));
 744
 745   for (i = 0; i < n; i++) hess[i] = hessdata + i * n;
 746   for (i = 0; i < k; i++) cgrad[i] = cgraddata + i * n;
 747
 748   choldecomp(hess, n, 0.0, &maxadd);
 749   /**** do something if not pos. definite ****/
 750
 751   for (i = 0, ldL = 0.0; i < n; i++) ldL += log(hess[i][i]);
 752
 753   if (k > 0) {
 754     /* forward solve for (L^-1) cgrad^T */
 755     for (l = 0; l < k; l++) {
 756       for (i = 0; i < n; i++) {
 757         if (hess[i][i] != 0.0) cgrad[l][i] /= hess[i][i];
 758         for (j = i + 1; j < n; j++) cgrad[l][j] -= hess[j][i] * cgrad[l][i];
 759       }
 760     }
 761
 762     /* compute sigma and stdev */
 763     for (i = 0; i < k; i++) {
 764       for (j = i; j < k; j++) {
 765         for (l = 0, hess[i][j] = 0.0; l < n; l++)
 766           hess[i][j] += cgrad[i][l] * cgrad[j][l];
 767         hess[j][i] = hess[i][j];
 768       }
 769     }
 770
 771     choldecomp(hess, k, 0.0, &maxadd);
 772     /**** do something if not pos. definite ****/
 773     for (i = 0, ldcv = 0.0; i < k; i++) ldcv += log(hess[i][i]);
 774   }
 775   else ldcv = 0.0;
 776
 777   *val = (n - k) * log(ROOT2PI) - ldL - ldcv;
 778   if (! *detonly) *val += f;
 779 }
 780
 781 #ifdef SBAYES
 782
 783 loglapfront(fvals, cvals, ipars, dpars, val)
 784      double *fvals, *cvals;
 785      MaxIPars *ipars;
 786      MaxDPars *dpars;
 787      double *val;
 788 {
 789   int detonly = FALSE;
 790
 791   loglapdet(fvals, cvals, ipars, dpars, val, &detonly);
 792 }
 793
 794 /************************************************************************/
 795 /**                                                                    **/
 796 /**                        First Order Moments                         **/
 797 /**                                                                    **/
 798 /************************************************************************/
 799
 800 moms1front(gf, n, m, x, hessdata, mean, stdev, sigmadata, h, typx)
 801      char *gf;
 802      int *n, *m;
 803      double *x, *hessdata, *mean, *stdev, *sigmadata, *h, *typx;
 804 {
 805   int i, j, k;
 806   double *hess, *sigma, *delg;
 807   double *delgdata, maxadd;
 808
 809   hess = (double **) S_alloc(*n, sizeof(double *));
 810   sigma = (double **) S_alloc(*m, sizeof(double *));
 811   delg = (double **) S_alloc(*m, sizeof(double *));
 812   delgdata = (double *) S_alloc(*m * *n, sizeof(double));
 813
 814   for (i = 0; i < *n; i++) hess[i] = hessdata + i * *n;
 815   for (i = 0; i < *m; i++) sigma[i] = sigmadata + i * *m;
 816   for (i = 0; i < *m; i++) delg[i] = delgdata + i * *n;
 817
 818   gevalfront(gf, n, m, x, h, typx, mean, delgdata);
 819
 820   /* get the cholesky decomposition L of the hessian */
 821   choldecomp(hess, *n, 0.0, &maxadd);
 822
 823   /* forward solve for (L^-1) delg^T */
 824   for (k = 0; k < *m; k++) {
 825     for (i = 0; i < *n; i++) {
 826       if (hess[i][i] != 0.0) delg[k][i] /= hess[i][i];
 827       for (j = i + 1; j < *n; j++) delg[k][j] -= hess[j][i] * delg[k][i];
 828     }
 829   }
 830
 831   /* compute sigma and stdev */
 832   for (i = 0; i < *m; i++) {
 833     for (j = i; j < *m; j++) {
 834       for (k = 0, sigma[i][j] = 0.0; k < *n; k++)
 835         sigma[i][j] += delg[i][k] * delg[j][k];
 836       sigma[j][i] = sigma[i][j];
 837     }
 838     stdev[i] = sqrt(sigma[i][i]);
 839   }
 840 }
 841
 842 /************************************************************************/
 843 /**                                                                    **/
 844 /**                        Second Order Moments                        **/
 845 /**                                                                    **/
 846 /************************************************************************/
 847
 848 typedef struct {
 849   MaxIPars max;
 850   int full, covar;
 851 } MomIPars;
 852
 853 typedef struct {
 854   MaxDPars max;
 855   double mgfdel;
 856 } MomDPars;
 857
 858 moms2front(ff, gf, gnum, x, typx, fvals, gvals, ipars, dpars,
 859            mean, stdev, sigmadata)
 860      char **ff, **gf;
 861      int *gnum;
 862      double *x, *typx, *fvals, *gvals, *mean, *stdev, *sigmadata;
 863      MomIPars *ipars;
 864      MomDPars *dpars;
 865 {
 866   char *msg;
 867   double h, loglap0, loglap1, loglap2;
 868   double *tilts, *fvals1, *gvals1, *x1;
 869   MomDPars dp1, *dpars1 = &dp1;
 870   MomIPars ip1, *ipars1 = &ip1;
 871   int i, n, m;
 872
 873   n = ipars->max.n;
 874   m = *gnum;
 875   h = dpars->max.h;
 876
 877   tilts = (double *) S_alloc(m, sizeof(double));
 878   x1 = (double *) S_alloc(n, sizeof(double));
 879   fvals1 = (double *) S_alloc(1 + n + n * n, sizeof(double));
 880   gvals1 = (double *) S_alloc(m + n * m, sizeof(double));
 881
 882   maxfront(ff, nil, nil, x, typx, fvals, gvals, nil, nil,
 883            ipars, dpars, nil, &msg);
 884   copypars(x, fvals, gvals, ipars, dpars, x1, fvals1, gvals1, ipars1, dpars1);
 885   loglapfront(fvals1, nil, ipars1, dpars1, &loglap0);
 886   copypars(x, fvals, gvals, ipars, dpars, x1, fvals1, gvals1, ipars1, dpars1);
 887   moms1front(gf, &n, &m, x1, fvals1 + n + 1, mean, stdev, sigmadata, &h, typx);
 888
 889   if (ipars->full) {
 890     for (i = 0; i < m; i++) {
 891       copypars(x, fvals, gvals, ipars, dpars,
 892                x1, fvals1, gvals1, ipars1, dpars1);
 893       ipars1->max.m = 1;
 894       ipars1->max.exptilt = FALSE;
 895       dpars1->max.newtilt = 1.0;
 896       maxfront(ff, gf + i, nil, x1, typx, fvals1, gvals1, nil, nil,
 897                ipars1, dpars1, nil, &msg);
 898       loglapfront(fvals1, nil, ipars1, dpars1, &loglap1);
 899       loglap1 -= loglap0;
 900
 901       copypars(x, fvals, gvals, ipars, dpars,
 902                x1, fvals1, gvals1, ipars1, dpars1);
 903       ipars1->max.m = 1;
 904       ipars1->max.exptilt = FALSE;
 905       dpars1->max.newtilt = 2.0;
 906       maxfront(ff, gf + i, nil, x1, typx, fvals1, gvals1, nil, nil,
 907                ipars1, dpars1, nil, &msg);
 908       loglapfront(fvals1, nil, ipars1, dpars1, &loglap2);
 909       loglap2 -= loglap0;
 910
 911       mean[i] = exp(loglap1);
 912       stdev[i] = sqrt(exp(loglap2) - exp(2.0 * loglap1));
 913       if (ipars->covar) sigmadata[i + i * m] = stdev[i] * stdev[i];
 914     }
 915     if (ipars->covar) {
 916       char *cgf[2];
 917       int j;
 918
 919       for (i = 0; i < m; i++) {
 920         for (j = i + 1; j < m; j++) {
 921           cgf[0] = gf[i];
 922           cgf[1] = gf[j];
 923           copypars(x, fvals, gvals, ipars, dpars,
 924                    x1, fvals1, gvals1, ipars1, dpars1);
 925           ipars1->max.m = 2;
 926           ipars1->max.exptilt = FALSE;
 927           dpars1->max.newtilt = 1.0;
 928           maxfront(ff, gf + i, nil, x1, typx, fvals1, gvals1, nil, nil,
 929                    ipars1, dpars1, nil, &msg);
 930           loglapfront(fvals1, nil, ipars1, dpars1, &loglap1);
 931           loglap1 -= loglap0;
 932
 933           sigmadata[i + j * m] = exp(loglap1) - mean[i] * mean[j];
 934           sigmadata[j + i * m] = sigmadata[i + j * m];
 935         }
 936       }
 937     }
 938   }
 939   else {
 940     for (i = 0; i < m; i++)
 941       tilts[i] = (stdev[i] > 0.0) ? dpars->mgfdel / stdev[i] : dpars->mgfdel;
 942
 943     for (i = 0; i < m; i++) {
 944       copypars(x, fvals, gvals, ipars, dpars,
 945                x1, fvals1, gvals1, ipars1, dpars1);
 946       ipars1->max.m = 1;
 947       ipars1->max.exptilt = TRUE;
 948       dpars1->max.newtilt = tilts[i];
 949       maxfront(ff, gf + i, nil, x1, typx, fvals1, gvals1, nil, nil,
 950                ipars1, dpars1, nil, &msg);
 951       loglapfront(fvals1, nil, ipars1, dpars1, &loglap1);
 952       loglap1 -= loglap0;
 953
 954       copypars(x, fvals, gvals, ipars, dpars,
 955                x1, fvals1, gvals1, ipars1, dpars1);
 956       ipars1->max.m = 1;
 957       ipars1->max.exptilt = TRUE;
 958       dpars1->max.newtilt = -tilts[i];
 959       maxfront(ff, gf + i, nil, x1, typx, fvals1, gvals1, nil, nil,
 960                ipars1, dpars1, nil, &msg);
 961       loglapfront(fvals1, nil, ipars1, dpars1, &loglap2);
 962       loglap2 -= loglap0;
 963
 964       mean[i] = (loglap1 - loglap2) / (2.0 * tilts[i]);
 965       stdev[i] = sqrt((loglap1 + loglap2) / (tilts[i] * tilts[i]));
 966       if (ipars->covar) sigmadata[i + i * m] = stdev[i] * stdev[i];
 967     }
 968     if (ipars->covar) {
 969       char *cgf[2];
 970       double ctilt, tscale[2];
 971       int j;
 972
 973       ctilt = dpars->mgfdel;
 974       for (i = 0; i < m; i++) {
 975         for (j = i + 1; j < m; j++) {
 976           cgf[0] = gf[i];
 977           cgf[1] = gf[j];
 978           tscale[0] = stdev[i] > 0 ? stdev[i] : 1.0;
 979           tscale[1] = stdev[j] > 0 ? stdev[j] : 1.0;
 980
 981           copypars(x, fvals, gvals, ipars, dpars,
 982                    x1, fvals1, gvals1, ipars1, dpars1);
 983           ipars1->max.m = 2;
 984           ipars1->max.exptilt = TRUE;
 985           dpars1->max.newtilt = ctilt;
 986           maxfront(ff, cgf, nil, x1, typx, fvals1, gvals1, nil, nil,
 987                    ipars1, dpars1, tscale, &msg);
 988           loglapfront(fvals1, nil, ipars1, dpars1, &loglap1);
 989           loglap1 -= loglap0;
 990
 991           copypars(x, fvals, gvals, ipars, dpars,
 992                    x1, fvals1, gvals1, ipars1, dpars1);
 993           ipars1->max.m = 2;
 994           ipars1->max.exptilt = TRUE;
 995           dpars1->max.newtilt = -ctilt;
 996           maxfront(ff, cgf, nil, x1, typx, fvals1, gvals1, nil, nil,
 997                    ipars1, dpars1, tscale, &msg);
 998           loglapfront(fvals1, nil, ipars1, dpars1, &loglap2);
 999           loglap2 -= loglap0;
1000
1001           sigmadata[i + j * m] = stdev[i] * stdev[j]
1002             * ((loglap2 + loglap1) /(2 * ctilt * ctilt) - 1.0);
1003           sigmadata[j + i * m] = sigmadata[i + j * m];
1004         }
1005       }
1006     }
1007   }
1008 }
1009
1010 static copypars(x, f, g, ip, dp, x1, f1, g1, ip1, dp1)
1011      double *x, *f, *g, *x1, *f1, *g1;
1012      MomIPars *ip, *ip1;
1013      MomDPars *dp, *dp1;
1014 {
1015   int i, n, m, nf, ng;
1016
1017   n = ip->max.n;
1018   m = ip->max.m;
1019   nf = 1 + n + n * n;
1020   ng = m + n * m;
1021
1022   for (i = 0; i < n; i++) x1[i] = x[i];
1023   for (i = 0; i < nf; i++) f1[i] = f[i];
1024   for (i = 0; i < ng; i++) g1[i] = g[i];
1025   *ip1 = *ip;
1026   *dp1 = *dp;
1027 }
1028
1029 /************************************************************************/
1030 /**                                                                    **/
1031 /**                          Laplace Margins                           **/
1032 /**                                                                    **/
1033 /************************************************************************/
1034
1035 lapmar1front(ff, gf, x, typx, fvals, ipars, dpars, xmar, ymar, nmar)
1036      char **ff, **gf;
1037      double *x, *typx, *fvals, *xmar, *ymar;
1038      MaxIPars *ipars;
1039      MaxDPars *dpars;
1040      int *nmar;
1041 {
1042   char *msg;
1043   int i, n, m, mindex;
1044   double h, loglap0, loglap1, xmode, stdev, sigmadata, ctarget[1];
1045   double *fvals1, *x1, *cvals, *cvals1, *fvals2, *x2, *cvals2;
1046   MaxDPars dp1, dp2, *dpars1 = &dp1, *dpars2 = &dp2;
1047   MaxIPars ip1, ip2, *ipars1 = &ip1, *ipars2 = &ip2;
1048
1049   n = ipars->n;
1050   m = 1;
1051   h = dpars->h;
1052
1053   x1 = (double *) S_alloc(n + 1, sizeof(double));
1054   fvals1 = (double *) S_alloc(1 + n + n * n, sizeof(double));
1055   cvals = (double *) S_alloc(m + n * m, sizeof(double));
1056   cvals1 = (double *) S_alloc(m + n * m, sizeof(double));
1057   x2 = (double *) S_alloc(n + 1, sizeof(double));
1058   fvals2 = (double *) S_alloc(1 + n + n * n, sizeof(double));
1059   cvals2 = (double *) S_alloc(m + n * m, sizeof(double));
1060
1061   maxfront(ff, nil, nil, x, typx, fvals, nil, nil, nil,
1062            ipars, dpars, nil, &msg);
1063   cpmarpars(x, fvals, cvals, ipars, dpars, x1, fvals1, cvals1, ipars1, dpars1);
1064   loglapfront(fvals1, nil, ipars1, dpars1, &loglap0);
1065   cpmarpars(x, fvals, cvals, ipars, dpars, x1, fvals1, cvals1, ipars1, dpars1);
1066   moms1front(gf, &n, &m, x1, fvals1 + n + 1, &xmode, &stdev, &sigmadata,
1067              &h, typx);
1068
1069   ipars->k = 1;
1070   ipars->vals_suppl = FALSE;
1071   ctarget[0] = xmode;
1072   maxfront(ff, nil, gf, x, typx, fvals, nil, cvals, ctarget,
1073            ipars, dpars, nil, &msg);
1074
1075   for (mindex = 0; mindex < *nmar && xmar[mindex] <= xmode; mindex++);
1076
1077   cpmarpars(x, fvals, cvals, ipars, dpars, x1, fvals1, cvals1, ipars1, dpars1);
1078   for (i = mindex; i >= 0; i--) {
1079     ctarget[0] = xmar[i];
1080     maxfront(ff, nil, gf, x1, typx, fvals1, nil, cvals1, ctarget,
1081              ipars1, dpars1, nil, &msg);
1082     cpmarpars(x1, fvals1, cvals1, ipars1, dpars1, x2,
1083               fvals2, cvals2, ipars2, dpars2);
1084     loglapfront(fvals2, cvals2, ipars2, dpars2, &loglap1);
1085     ymar[i] = exp(loglap1 - loglap0);
1086   }
1087   cpmarpars(x, fvals, cvals, ipars, dpars, x1, fvals1, cvals1, ipars1, dpars1);
1088   for (i = mindex + 1; i < *nmar; i++) {
1089     ctarget[0] = xmar[i];
1090     maxfront(ff, nil, gf, x1, typx, fvals1, nil, cvals1, ctarget,
1091              ipars1, dpars1, nil, &msg);
1092     cpmarpars(x1, fvals1, cvals1, ipars1, dpars1, x2,
1093               fvals2, cvals2, ipars2, dpars2);
1094     loglapfront(fvals2, cvals2, ipars2, dpars2, &loglap1);
1095     ymar[i] = exp(loglap1 - loglap0);
1096   }
1097 }
1098
1099 static cpmarpars(x, f, g, ip, dp, x1, f1, g1, ip1, dp1)
1100      double *x, *f, *g, *x1, *f1, *g1;
1101      MaxIPars *ip, *ip1;
1102      MaxDPars *dp, *dp1;
1103 {
1104   int i, n, k, nf, ng;
1105
1106   n = ip->n;
1107   k = ip->k;
1108   nf = 1 + n + n * n;
1109   ng = k + n * k;
1110
1111   for (i = 0; i < n; i++) x1[i] = x[i];
1112   for (i = 0; i < nf; i++) f1[i] = f[i];
1113   for (i = 0; i < ng; i++) g1[i] = g[i];
1114   *ip1 = *ip;
1115   *dp1 = *dp;
1116 }
1117 #endif /* SBAYES */
1118
1119 /*
1120   TODO
1121
1122   get hessian from gradiant for analytical gradiants
1123
1124   avoid repeated derivative calls in mimimize.
1125
1126   2d margins
1127
1128   use pos. definiteness info in margins
1129
1130 */